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Abstract 

pde2path is a free and easy to use Matlab continuation/bifurcation package for elliptic systems 
of PDEs with arbitrary raany components, on general two dimensional domains, and with rather 



Rayleigh-Benard convection, and von Karman plate equations. These serve as templates to 
Oh study new problems, for which the user has to provide, via Matlab function files, a descrip- 

^ tion of the geometry, the boundary conditions, the coefficients of the PDE, and a rough initial 

guess of a solution. The basic algorithm is a one parameter arclength-continuation with op- 
tional bifurcation detection and branch-switching. Stability calculations, error control and 
mesh-handling, and some elementary time-integration for the associated parabolic problem are 
,_, also supported. The continuation, branch-switching, plotting etc are performed via Matlab 

pLn command-line function calls guided by the AUTO style. The software can be downloaded from 

www.staff.uni-oldenburg.de/hannes.uecker/pde2path, where also an online documenta- 
^ tion of the software is provided such that in this paper we focus more on the mathematics and 

the example systems. 
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1 Introduction 

For algebraic systems, ordinary differential equations (ODEs), and partial differential equations 
(PDEs) in one spatial dimension there is a variety of software tools for the numerical continu- 
ation of families of equilibria and detection and following of bifurcations. These include, e.g., 
AUTO [8], XPPaut [7] (which rehes on AUTO for the continuation part) and MatCont [11], see also 
www.enm.bris.ac.uk/staff/hinke/dss/ for a comprehensive though somewhat dated list. An- 
other interesting approach is the "general continuation core" coco, [28]. 

However, for elliptic systems of PDEs with two spatial dimensions there appear to be few general 
continuation/bifurcation tools and hardly any that work out-of-the-box for non-expert users. ^ Our 
software pde2path is intended to fill this gap. Its main design goals and features are: 

• Flexibility and versality. The software treats PDE systems 

G{u, A) := - V • (c (g) Vu) + au - b ®Vu - f = 0, (1) 

where u = u{x) G M^, x ^ ft C M? some bounded domain, A G M is a parameter, c G 
rNxNx2x2^ 6 G R^x^x2 (see (4), (5) below), a G and / G can depend on x, u, Vu, 

and, of course, parameters.^ The current version supports "generalized Neumann" boundary 
conditions (BC) of the form 

n • (c (g) Vu) + qu = g, (2) 

^PLTMG [2] treats scalar equations, and there are many case studies using ad hoc codes, often based 
on AUTO using suitable expansions for the second spatial direction; for 2D systems there also is ENTWIFE, 
www. sercoassurance . com/entwife/introduction. html, which however appears to be no longer maintained since 
2001. For experts we also mention Loca [26], which is designed for large scale problems, and oomph [13], another large 
package which also supports continuation/bifurcation, though this is not yet documented. 

^The standard assumption is that c, a, /, 6 depend on u,Vu, . . . locally, e.g., f(x,u) = f{x,u{x)); however, the 
dependence of c, a, /, b on arguments can in fact be quite general, for instance involving global coupling, see §3.5. In 
particular, we added the —b Vu term to the pdetoolbox-form for the effective evaluation of Jacobians, see below. 
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where n is the outer normal and again q G R and g G M may depend on Vu 
and parameters. These boundary conditions include zero flux BC, and large prefactors in 
g generate a "stiff spring" approximation of Dirichlet BC that we found to work reasonably 
well. 

There are a number of predefined functions to specify domains and boundary conditions, 
or these can be exported from matlab's pdetoolbox GUI, thus making it easy to deal with 
(almost) arbitrary geometry and boundary conditions. 

The software can also be used to time-integrate parabolic problems of the form 

dtu = -G{u,X), (3) 

with G as in (1). This is mainly intended to easily find initial conditions for continuation. 
Finally, any number of eigenvalues of the Jacobian Gu{u^X) can be computed, thus allowing 
stability inspection for stationary solutions of (3). 

• Easy usage. The user has to provide a description of the geometry, the boundary conditions, 
the coefficients of the PDE, and a rough initial guess of a solution. There are a number of 
templates for each of these steps which cover some standard cases and should be easy to adapt. 
The software provides a number of Matlab functions which are called from the command line 
to perform continuation runs with bifurcation detection, branch switching, time integration, 
etc. 

• Easy hackability and customization. While pde2path works "out-of-the-box" for a 
significant number of examples, already for algebraic equations and ID boundary value prob- 
lems it is clear that there cannot be a general purpose "solve-it-all" tool for parametrized 
problems, see, e.g., [30, Chapter 3]. Thus, given a particular problem the user might want to 
customize pde2path. We tried to make the data structures and code as modular and trans- 
parent as possible. When dealing with a tradeoff between speed and readability we usually 
opted for the latter, and thus we believe that the software can be easily modified to add new 
features. In fact, we give some examples of "customization" below. Here, of course, having 
the powerful Matlab machinery at our disposal is a great advantage. 

Remark 1.1. The i^^ components of V • (c V?i), au and b ® Vu in (1) are given by 

N 

[V • (c Vu)]i := ^[dxCijiidx + d^Cijudy + dyCij2idx + dyCij22dy]uj , (4) 

N N 

[au]i = ^QijUj, [b (g) Vu]i := ^[hjidx + bij2dy\uj, (5) 

and / = (/i, . . . , /tv) should be seen as a column vector. If, for instance, we want to implement 
—DAu = —{diAui^ . . . , (InAun) = —V • {D\/u) with D a constant diagonal diffusion matrix, as it 
often occurs in applications, then 

Cull = Cii22 = di, i = 1, . . . , TV, and all other Cijki = 0, (6) 

and there are special ways to encode this (and other symmetric situations for c and a) in the 
pdetoolbox. See the templates below, and §3.1.3. For c^a^f see also the pdetoolbox documenta- 
tion, for instance assempde in the Matlab help. pde2path also provides a simplified encoding for 
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isotropic systems without mixed derivatives, see §4.1. Finally, the i*^ component of n • (c® Vu) is 
given by 

N 

[n • (c (g) Vu)]i = "^[niicijiidx + Cijudy) + n2{cij2idx + Cij22dy)]uj, (7) 

where n = (^1,712), J 

Remark 1.2. Clearly, the splitting between a and / (6 and /) in (3) is not unique, e.g., for 
G{u) = —Au — Xu + we could use (a — —A, / — —u^) or (a = 0, / = A?i — Similarly, for, 
e.g., G{u) = — A^i — 9^^!^ we can use b = (1, 0) and / = or & = (0, 0) and / = dxU. This flexibility of 
(1) has the advantage that in most cases the needed derivatives G^, Gx can be assembled efficiently 
from suitable coefficients c, a, 6, and no numerical Jacobians are needed. 

Also note that (1) allows to treat equations in nondivergence form, too. For instance, we may 
write a scalar equation —c{u)Au — f{u) = as —V • {c{u)Vu) + {c'{u)Vu) • Vu — f{u) = 0, and set 
biii{u) = —c'{u)dxU and &ii2('i^) = —c'{u)dyU^ or add —{c'{u)Vu) • Vu to /. J 

Currently, the main drawbacks of pde2path are: 

• pde2path requires Matlab including the pdetoolbox. Its usage explains the form (1). One 
of its drawbacks is a somewhat slow performance, compared to, e.g., some Fortran imple- 
mentations of the FEM.^ On the other hand, in addition to the Mat lab-environment, the 
pdetoolbox has a number of nice features: it also takes care of the geometry and mesh gen- 
eration, it is well documented, it is fully based on sparse linear algebra techniques (which are 
vital for large scale problems), it exports (sparse) mass and stiffness matrices, and it provides 
a number of auxiliary functions such as adaptive mesh-refinement, or various plot options. 

• Presently, only one parameter continuation is supported, and only bifurcations via simple 
eigenvalues are detected, located, and dealt with^. We plan to add new features as examples 
require them, and invite every user to do so as well. 

In the following we first very briefiy recall some basics of continuation and bifurcation. Then 
we explain design and usage of our software by a number of examples, mainly a modified Bratu 
problem as a standard scalar elliptic equation, some AUen-Cahn type equations, some pattern 
forming Reaction-Diffusion systems, including some animal coats intended for illustration of how 
to set up problems with complicated geometries. We give a rather detailed bifurcation diagram for 
the Schnakenberg system, and we consider three rather classical problems from physics: Rayleigh- 
Benard convection, some multi-component Bose-Einstein systems, and the von Karman plate 
equations. Thus, besides some mathematical aspects of continuation and the example systems, 
here we explain the syntax and usage of the software in a rather concise way. More comprehensive 
documentation of the data structures and functions is included in the software, or online at [23]. 

Acknowledgements. We thank Uwe Priifert for providing his extension of the pdetoolbox and 
the documentation. Users familiar with AUTO will recognize that AUTO has been our guide in 
many respects, in particular concerning the design of the user interface. We owe a lot to that 
great software. We also thank Tomas Dohnal for testing early versions of pde2path and providing 
valuable hints for making pde2path and this manual more user friendly. JR acknowledges support 
by the NDNS+ cluster of the Dutch Science Fund NWO. 

^Another drawback is a somewhat unhandy non-GUI description of geometry and boundary conditions, but for 
these we provide fixes. See also, e.g., [24]. 

^In case symmetries cause multiple eigenvalues, artificial symmetry breaking sometimes is a viable ad hoc solution 
for the latter 
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Remark 1.3. Please report any bugs to pde2path@mathematik.uni-oldenburg.de, as well de- 
sired additional features. We will appreciate any feedback, and will be happy to provide help with 
interesting applications. See also [23] for an online documentation of the software, updates, FAQ, 
and general further information. 



2 Some basics of continuation and bifurcation 

2.1 Arclength continuation 

A standard method for numerical calculation of solution branches of A) = 0, where G : 
X X M X is at least C""^, X a Banach space, is (pseudo) arclength continuation. For convenience 
and reference here we recall the basic ideas. Standard textbooks on continuation and bifurcation 
are [30, 12, 17, 1], see also [16, 5], and the "matrix-free" approach [9]. Consider a branch z{s) := 
{u{s)^ X{s)) G X X M parametrized by 5 G M and the extended system 

H{u, \)=( = G X X R, (8) 



p{u,\s] 

where p is used to make s an approximation to arclength on the solution arc. Assuming that 
X is a Hilbert space with inner product (•,•)? the standard choice is as follows: given and 
a point (?io5^o) '-— ('^(•^o)? ^(-^o))? ^.nd additionally knowing a tangent vector tq := (^io^-^o) •= 
^^{u{s),X{s))\s=so we use, for s near sq, 

p{u, A, s) := i {iio, u{s) - uq) + (1 - O^oiK^) - Aq) - {s - sq). (9) 

Here < ^ < 1 is a weight, and tq is assumed to be normalized in the weighted norm 

\\r\k ■■= { (a) ' Q )^ -■= ^ ^) + (1 - 

For fixed s and ||to||^ = 1, p{u^X^s) = thus defines a hyperplane perpendicular (in the inner 
product (•, •)^) to To at distance ds := 5 — sq from (uq, Aq). We may then use a predictor (ix^, A^) = 
{uq^ Aq) +ds To for a solution (8) on that hyperplane, followed by a corrector using Newton's method 
in the form 

= © ""^'^-'KS (l-«)Ao)- 

Since dgP = — 1, on a smooth solution arc we have 

Thus, after convergence of (10) yields a new point (?ii, Ai) with Jacobian A^, the tangent direction 
Ti at (ui^Xi) with conserved orientation, i.e., (ro,ri) = 1, can be computed from 



A^Ti = , with normalization ||ti||^ = 1. 



(12) 



Alternatively to (10) we may also use a chord method, where A = A{u^^X^) is kept fixed during 
iteration, 

(a'+i) = (a') (13) 
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This avoids the costly evaluation of Gu at the price of a usually modest increase of required 
iterations. 

The role of ^ is twofold.^ First, if G(u, A) = comes from the discretization of a PDE G{u, A) 
such as (1) over a domain Q with Up spatial points, then u G with large say p = Nrip. 



Typically, we want to choose ^ such that ^||u||^p is an approximation of 
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If (as usual). 



P|ll"llL2(n)- 

u = 1 corresponds to Uj = 1 for j = 1, . . . , n^, then a rough estimate can be obtained by assuming 
that each component = 1, z = 1, . . . , A^. Then 



1 

W\ 



\u\ 



N = CII^IIrp — hence ^ = l/rip 



(14) 



This gives the basic formula for our choice of ^. It is important that different ^ may give different 
continuations: in the Newton loop, small ^ favors changes in while larger ^ favors A, see Fig.l 
for a sketch. Moreover, ^ is also related to the scaling of the problem: if, e.g., we replace A by, 
say, A := 100 A, then ^ should be adapted accordingly, i.e., ^ = summary, ^ should be 

considered as a parameter that can be used to tune the continuation, and that may also be changed 
during runs if appropriate. 







5=0.1 






y 5=0.9 






l< Z(S) (3) 

\ 1 






\ 5=0.5 



Figure 1: The role of ^ in a one-dimensional (u G M) sketch with tq = ('Uo^Aq) = (1,1) (unnormalized). 
Depending on ^ we get different hyperplanes {(ii,A) G : (ro,(ii,A))^ = ds} and consequently different 



next points 



1,2,3, on the solution curve z{s). Small ^ favors Newton search in u direction (and 



thus orthogonal to a "horizontal" branch), while large ^ favors the A direction (parallel to a "horizontal" 
branch) . 

Given a weight ^, a starting point (uq^ Aq, tq), and an intended step size ds, the basic continuation 
algorithm thus reads as follows, already including some elementary stepsize control: 



Algorithm cent 

1. Predictor. Set {u^,X^) = (uq^Xq) + dsTo. 

2. Newton— corrector. Iterate (10) (or (13)) until convergence; decrease ds if (10) fails to 
converge and return to 1; increase ds for the next step if (10) converges quickly; 

3. New tangent. Calculate ri from (12), set {uo^Xq^tq) = (t^i, Ai,ti) and return to step 1. 



Theoretically, this does not work at possible "bifurcation points" where A is singular.^ More 
specifically, we define: 



^Here u stands for the FEM approximation of and G(u,A) for the FEM approximation of G. Below, the 
difference between the two will be clear from the context, and thus we will mostly drop the different notations again. 
^Allthough generically continuation routines simply shoot past singular points. 



6 



Bl. A simple bifurcation point is a point {u^ A) where detA changes sign. The imphcit assumption 
is that this happens due to a simple eigenvalue of A crossing zero. 

This clearly excludes folds (aka turning points), where a simple eigenvalue of A reaches zero 
but detA does not change sign, see [16]. However, folds are no problem for the algorithm and can 
easily seen in the bifurcation diagram anyway. Therefore there is no special treatment of folds in 
the current 1-parameter version of pde2path. Bl also excludes bifurcations via even numbers of 
eigenvalues crossing, which are more complicated to deal with. 

Remark 2.1. Numerically, for Bl we found it more robust to use ^ = 1/2 in the definition of A for 
bifurcation purposes. For algorithmic reasons, we only use the first part of Bl for the detection of 
bifurcation points, i.e., the sign change of detA, which also occurs for an odd number of eigenvalues 
(counting multiphcities) crossing zero. Finally, to calculate sign(det A) we only calculate the set 
So of eigenvalues /i^, i = 1, . . . , ngig of A closest to and then use 

sign( det A)=sign(n^5Re/i^). (15) 

Here the implicit assumption is that ngig is sufficiently large such that all eigenvalues with negative 
real part are always contained in Sq. This is reasonable as + 7 is a positive elliptic operator for 
sufficiently large 7.^ J 

After detection of a bifurcation between Sk and s^^i, the bifurcation is located by a bisection 
method. To switch branches we use "Method I" of [16] (page 379). Let (uq^Xo) be a simple 
bifurcation point, Gu = ^^('^05-^0)5 ^.nd tq = (uq^Xq) be the tangent along the branch already 
computed. To obtain a tangent ri along the other branch we proceed as follows: 

Algorithm swibra 

1. Calculate 01,^1 with G^^i = 0,G^^^i = 0, = 1, {^1,^1) = 1. 

2. Let ao = Aq, ai = ip^uo, 0o = c^o^iuo - <^i0i). 

3. Choose some small 5 > and calculate the finite differences 

ai = ^V^f [G^(n + 50i,Ao) -G^^i, 

[Gu{uo + 6^i,Xo) - G^]^o + Gx{uo + 6^1, Xo) -Ga(^o,Ao) . 

Assuming o^o 7^ (see [16] if this is not the case), set 

a, = - + 26i) and n = ("^'^\+ "^'^°) • 

Choose^ a weight ^ and a stepsize ds, set tq = ti/||ti||^ and go to cont, step 1. 
^ If branch-switching fails, i.e., if there is no convergence in cont or if the solution falls back onto the known 
branch, then it may help to change ds and/or ^. 

^Due to occasional numerical problems in the eigenvalue calculations, in the current standard setting of pde2path 
(controlled by switches, see below) we actually combine the sign-change of det (A) with a consistency check with the 
eigenvalues of Gu- Moreover, in practice it is sometimes useful to turn off branch point detection and localization in 
the initial continuation (again by setting switches, see below), but only monitor changes in the number of eigenvalues 
with negative real part in Gu- For each such change select a nearby starting point for a new continuation with branch 
detection. See also §5.2. Finally, we also provide a routine f indbif which first scans a branch for a change of the 
number of unstable eigenvalues, and then uses Bl to locate a bifurcation. 

Our standard setting is rieig = 50, but of course this is highly problem dependent and should be adapted by the 
user when needed. We give a warning if > |/Xneig|/2 since then eigenvalues might wander out of Eo to the left in 
the next steps. 



1 rr. 
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2.2 Switching back and forth to the natural parametrization 

If A := dgX does not change sign, then we know that a branch also has the "natural parametrization" 
(i/(A),A), and, except at possible bifurcation points, Gu{u^X)u' {X) = —G\{u^X) has the unique 
solution 

u'{\) = -Gu{u,\)-^Gx{u,\). 

Thus, given {u{)^ Aq) we may use the predictor {v}^ Ai) = (?i05 ^o) + ds (?i^(Ao), 1) and then correct 
with fixed A = Ai. Algorithmically, however, we choose to keep the predictor [u^^ Ai) = {uq^ Aq) + 
dsTo, i.e., altogether, 

{u\ Ai) = (no, Ao) + dsro, u^^^ = - Gu{u\ \i)-^G{u\ Ai). (16) 

After convergence to (i^i,Ai) we calculate the new tangent via (12), and Bl can again be used 
as a check for bifurcation. Moreover, with tolj^ > a given tolerance, say tolj^ = 0.5, this gives a 
criterion when to switch back and forth between the algorithms, namely: 

If |A| > tol^^, then use (16), else use (10). (17) 

Here again the weight ^ is important: for fixed ^ = 1/2 (say), A ^ as ^ oo unless the branch 
is strictly horizontal, i.e., = 0. 

Remark 2.2. If applicable, (16) is usually slightly faster than (10), as expected. On the other 
hand, we found that even for "nearly horizontal" branches, locating the bifurcation point typically 
works better with arclength continuation (10). J 

3 Some scalar problems in pde2path 

We now start the tutorial on pde2path by way of basic examples. The names in brackets refer to 
the sub-directory name of the directory demos, which contains the given example. 

3.1 Bratu's problem (bratu) 

Our first example is the scalar elliptic equation 

-/\u - f{u, A) = 0, f{u, A) = -l{){u - Ae^), u = u{x) e M, (18) 
on the unit square with zero flux BC, i.e., 

xen = (-1/2, 1/2)2, ^^^1^^ _ Q ^19) 

This problem has the advantage that a number of results can immediately be obtained analytically, 
that there are some nontrivial numerical questions (see below), and that we can compare with 
previous results, see, e.g., [2, 3, 21] 

There is a primary homogeneous branch u = Uh{s)^ X = X{s)^ "starting" in (0,0), on which 
{uh{s)^X{s)) satisfies f{u) = —10{u — Ae^) = 0. Bifurcation points {u^^Xk) are obtained from 

GuW —Aw — fuW = which yields 10{uk — 1) = /J^k where jj.^ {k\ + k^Tv'^, k G Nq, are the 
eigenvalues of —A on see Table 1. From §2.1, for arclength continuation the fold is nothing 
special and the two dimensional kernels k = (^1,^2), ki ^ A;2, will go undetected using Bl. The 
simple bifurcation points should be detected, and branch switching can be tried. On bifurcating 
branches, further bifurcations may be expected. 

In pde2path, (18), (19) can be setup and run in a few steps explained next, to quickly obtain 
the (basic) bifurcation diagram and solution plots in Fig. 2 on page 12. 
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k 


(0,0) 


(1,0), (0,1), 


(1,1) 


(2,1),(1,2) 


(2,2) 




Uk 

type 


1 

1/e « 0.3679 
fold 


l + vrViO 
« 0.2724 

double 


1 + 7rV5 
« 0.1520 
simple 


1 + Try 2 
^ 0.0157 
double 


1 + 47rV5 
« 0.0012 
simple 





Table 1: Bifurcation from homogeneous branch in (18). 



3.1.1 Installation and preparation 

The basic pde2path installation consists of a root directory, called pde2path, with a subdirectory 
p2plib containing the actual software, a subdirectory demos with a further subdirectory for each 
problem, a subdirectory octcomp, providing some basic octave compatibility (see [23]), and one 
Mat lab file setpde2path.m, which is a utility function to set the Mat lab path. Each of the demos 
comes with a file *cmds .m, which contains the commands to run the example (and some comments), 
and which should be seen as a quarry for typical commands, and with a file *demo .m, which produces 
more verbose output. To start we recommend to run setpde2path (without arguments) in the root 
directory pde2path and then change into one of the demo-directories, e.g., type cd demos/bratu in 
Matlab. Then inspect the file bratucmds . m and copy paste the commands to the Matlab command 
line, or just execute bratucmds or bratudemo. 

3.1.2 General structure, initialization, and continuation runs 

In pde2path, a continuation and bifurcation problem is described by a structure, henceforth called 
p (as in problem), which we now outline, see also Tables 2 and 4.^ Essentially, p contains 

• function handles which describe the functions c, a, 6, / and the BC (and possibly the Jacobian) 
in (1); 

• fields which describe the geometry of the problem, including the FEM mesh; 

• fields which hold the current solution, i.e., A and the tangent r; 

• a number of variables, switches and further functions (i.e., function handles) controlling the 
behaviour of the continuation and bifurcation algorithm, and filenames for file output. 

Studying a continuation and bifurcation problem using pde2path thus consist of: 

• Setting up a file defining the coefficient functions c, a, 6 and / (and usually a function for the 
Jacobians) in (1), e.g. bratuf .m (and bratujac.m). (Here we assume that the BC function 
p.bcf is defined inline as in bratuinit .m and many of the further examples) 

• Setting up an initialization function file, e.g., bratuinit filling p. The main steps are (1) 
define p.f and p.jac, (2) define geometry and mesh, and usually the BC by an inline function, 
(3) set the parameters and provide a starting point. In many cases most parameters and 
switches can be set to "standard values". For this we provide the function p=stanparam(p) , 

^In addition to the fields/variables listed, there are quite a few more within p. See stanparajn.m for these, and 
also for more comments on the ones which are listed. Some of the "control fields" are unlikely to be changed by 
the user, at least at the beginning, e.g., p. evopts .disp=0 (to suppress output during eigenvalue calculations), or 
p . pf ig= 1 , p . brf ig=2 (the figure numbers for plotting) , and some of the additional fields / variables are only generated 
during computation, e.g., the residual res and the error estimate err, which are put into p for easy passing between 
subroutines and user access. Currently there are no global variables in pde2path, with the exceptions pj ,lainj which 
are set for numerical differentiation in resinj, and possibly LU preconditioners for iterative linear system solvers, 
see §3.1.6. See also §3.6. 
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func. handles 


meaning 


example (in initbratu) 


f 

jac 

bcf 
outfu 

ufu 

headfu 
blss, Iss 


[c,a,f ,b]=f (p,u,lain) PDE coefficients in (1) 
[c,fu,f lam,b] =jac(p,u,lam) used to build Gu,Gx from 
/n, fx if desired (see jsw below), 
bc=bcf (p,u,lain), boundary conditions function 
out=outfu(p,u,lam) , defining output for bifurcation dia- 
gram, i.e., quantities saved on p.brcLnch. Default setting 
out=stanbra(p,u,lam) gives out=[||i//i ||oo5 ll'^illz/^]^ 
cstop=ufu(p,brout,ds), user function, called after each cal- 
culation of a point, e.g. for printing information and checking 
stopping criteria 

headfu (p) (no return arguments); print headline of screen out- 
put 

(bordered) linear system solver, see §3.1.6 


p.f=@bratuf 
p.jac=@bratujac 

p.bcf=@(p,u,lam) bc"'^ 
p.outfu=@stanbra^ 

p.ufu=@stanufu^ 

p.headfu=@stanheadfu^ 
p.blss=@blss^ p.lss=@lss^ 


various var. 


meaning 


example (in initbratu) 


geo 

points, edges, tria 
bpoints,..,btria 
neq, np, nt 
u,lam,tau,ds 
branch, bifvals 


geometry (and also BC in call to recnbcl) 
mesh 

background mesh (used for mesh- adaption) 
number of equations, mesh points, triangles 
essential continuation data 

lists generated via p. outfu, used to plot bif. diagrams 


[p.geo,bc]=recnbcl(0.5,0.5);"'^' 

llr O 7 J V ' / ' 

p=stanmesh(p,30,30) 
p=setbmesh(p) ^ 
automatically from mesh 
p.u=0.2*ones(p.np); . . . ; 
p.branch=[]; p.bifvals=[]; 


file names 


meaning 


standard (in setfn) 


pre 

pname,bpname 


name of subdir for files; by default automatically set to the 
struct name used in call to cont or *init 
(base) filenames for output of points, bifurcation points; actual 
filenames augmented by counter 


p.pname='p' for |)oint, 
p.bpname='bp' for 6if. point 



" typical setup with be independent of A, u and hence defined in advance 
^ "standard" choices provided by pde2path 

^ see bratuinit.m and documentation of pdetoolbox for explanation 

Table 2: Basic variables in structure p of a problem, and their initialization in bratuinit. 



main functions 


purpose 


p=cont(p) 

p=swibra(pre,ffie,npre) 

p=meshadac (p , varargin) 

p=findbif(p,ichange) 

p=loadp(pre,file,npre) 

plotbra(p,wnr,cmp,aux) 

plotbraf(pre,file,wnr,cmp,aux) 

plotsol(p,wnr,cmp,pstyle) 


main continuation routine 

branch-switching at bifurcation point from previous run, prefix set to npre 
adapt mesh in p, see §3.1.5 for varargin 

locate bifurcation point based on Gu'i cont itself uses bifdetec (bifchecksw> 0) 
load solution struct p from file and reset prefix to npre 
plot component cmp of branch over A, in figure number wnr 
as plotbra but from ffie (saved previously, leave out the '.mat'), 
plot solution, use plotsolf(pre,file,wnr,. . . ) to plot from file 



Table 3: Main "user" functions in pde2path. 



which should be called first, and afterwards individual parameters can be reset as needed. For 
the mesh generation we provide an elementary function p=stanmesh(p,hmax) , where hmax 
is the maximal triangle side-length. This is based on initmesh from the pdetoolbox, i.e., 
a Delaunay algorithm. For rectangular domains the syntax p=stanmesh(p,nx,ny) is also 
allowed, which is based on poimesh from the pdetoolbox, with obvious meaning. ^ 

• Calling a number of pde2path functions. The basic call is p=cont(p); (a continuation 
run), which can be followed, e.g., by a repeated call to extend the branch. Or, in case 
a branch point has been found, a call to branch switching and subsequent continuation 

^If applicable, poimesh obviously is faster and gives more regular meshes, if used with care, i.e., choose nx/ny 
according to Lx/Ly where Lx,Ly are the sidelengths of the rectangle. However, we found that in some cases the 
Delaunay mesh gives more robust numerics, in particular after mesh refinement (see §3.1.5). Thus we recommend to 
experiment with both. 
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Newton&Cont 


meaning 


imax,normsw 

tol 

nsw 

jsw 

dsmin , dsmax , dlammax 

icllllllllil , iclllllllolX 

nsteps 

parasw,lamdtol 
amod, maxt, ngen 
errchecksw, errtol 


Newton controls: max number of steps and selection of a norm ('norm-switch') 
stop-crit. for Newton, typically should be around le-10; 
for Newton, 1 for chord 

switch for derivatives (G^, Gx). 0: (Gu, Gx) by (c, fu,b, fx), 
1: Gu by c, A, b, Gx by FD, 2: Gu by FD, Gx by /a, 3: both by FD. 
min/max stepsize in s, max stepsize in A 

111111/ iiicLx A, piebeu uo -|--LU , leoet to uoe cib o topping ciiLeiici 
number of steps to take 

parametrization switch and tolerance: if parasw=0 resp. parasw=2 then always use 
(16) resp. (10). If parasw=l then use (17) with tolj^ =lamdtol 

controls for mesh adaption: adapt every amod-th step, aim at maxt triangles, in at 
most ngen refinement steps 

switch and tol for a posteriori error estimate and handling, see §3.1.5. 


Bif . ,Plot&User-control 


meaning 


bifchecksw 
neig 

eigsstart 
bisecmax 
spcalcsw 
pmod/pstyle 
pciiip, upciiip 

isw,vsw 
pfig,brfig,ifig 

timesw 
nbp 


for no checks, 1 for check via Bl with consistency with eigenvalues of Gu, 
2 for Bl alone 

number of eigenvalues to be calculated in spcalc and bif detec, default 50 
to start eigs randomly, 1 to start with (1, . . . , 1) (default) 
max number of bisections to locate bifurcation; turned off by bisecmax=0 
0/1 for eigenvalue calculation off /on 

plot each pmod-th step in style pstyle (1 mesh, 2 pcolor, 3 rendered 3D, . . . ) 
cuiiipuiieiiL Lu piuu, cuiiipoiieiii oi uicLiiciieo lo piou 

DcLVc; CVClj' oliUJH Lll oLCJJ 

interaction/verbosity switch: O=none, l=some, 2=much; 

figure- numbers for u-plot, branch-plot and info-plot during cont.; in ifig we plot add. in- 
formation, e.g., after mesh adaption, or the new tangent after swibra 
if > 0, print timing info after cont. See stanparam.m for the timers in p. 
number of user-components of branch to be printed on screen (p.ufu=@stanufu) 


resfac, mst, pmimax 


pmcont only, see §4.3 



Table 4: Main switches and controls in a structure p used in cont, see stanparam.m for typical 
values and a number of additional switches with detailed comments. See also §4.3 for additional 
parameters controlling pmcont, the parallelmulti-continuation version. 

by q=swibra(^p\ 'bpl\ ^qO ; q=cont(q); where, e.g., ./p/bpl.mat is data written at 
a branch point during the previous run. Inbetween runs, p can be modified from the com- 
mand line, e.g., type p.imax=5 (say) to (re)set the maximal number of Newton-iterations, or 
call p=meshref (p) to refine the mesh before a subsequent run; afterwards, call p=cont(p) 
again. According to the settings, data is plotted and written to files in a sub-directory with 
name p. pre. There are also functions for further postprocessing, e.g., plotting of solutions 
and bifurcation diagrams, whose documentation is mainly provided within the corresponding 
matlab files, and by the calls in the example directories. 

The fundamental user provided functions thus are the coefficient function p . f and the addition- 
ally recommended jacobian function p. jac, see §3.1.3. To study a new problem, we recommend to 
edit copies of the files *init.m, *f .m and *jac.m of a suitable example (e.g. *=bratu for a scalar 
problem) in an empty directory, and start with calling p=[] ;p=newinit (p) ; p=cont(p). The bi- 
furcation diagram in Fig. 2 and the solution plots are generated from the commands in Table 5, 
either given from the command line, or put into a Matlab script 

^°The figures generated from the commands in bratucmds.m (and similarly for the remaining demos to come) 
may differ slightly from the ones in this manual; this is due to minimal postprocessing via Matlab's "Edit Figure 
Properties" to set fontsizes, axis labels, etc. 
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function p=bratuinit (p) % init-routine , see bratuinit.m for more comments 
p=stanparam(p) ; p.neq=l; p.f=@bratuf; p. jac=@bratujac; [p . geo ,bc] =recnbcl (0 . 5 , . 5) ; 
p.bcf =@(p,u,lam) be; % typical inline definition of the BC function 
nx=20;p=stanmesh(p,nx,nx) ; p=setbmesh(p) ; % mesh and "background" mesh 
pre=sprintf ( ^7oS ^ , inputname(l) ) ; p=setf n(p,pre) ; % set filename (prefix) 
p.xi=l/p.np; p . dlammax=0 . 02 ; p . lammin=0 . 02 ;p . tau=l ; % set p.tau to something 
p.lam=0.2; p.u=0. l*ones(p.np, 1) ; p.ds=0.05; % "trivial" branch 

function [c,a,f ,b] =bratuf (p,u,lam) coeff for Bratu 
u=pdeintrp(p.points ,p. tria,u) ; c=l; a=0; f =-10* (u-lam*exp(u) ) ; b=0; 

function [c,a,f ,b]=bratujac(p,u,lam) %% Jacobian for Bratu 

u=pdeintrp(p.points ,p. tria,u) ; c=l; f u=-10* (l-lam*exp(u) ) ; f lam=10*exp(u) ; b=0; 

7o commands to run bratu in pde2path (selection, see also script file bratucmds.m) 
p= [] ; p=bratuinit (p) ; p=cont (p) ; 

q=swibra('p' , ^bpl^ , ^q') ; q.lammin=0. 1 ;q.nsteps=20;q=cont (q) ; 
plotbra(p,3,2, ^msS 12, ^lwS5, 'f s' , 16, 'cl\ 'k' ) ;plotsolf ( ^q^ , 'p20^ ,4,1,1) ; 



Table 5: The basic init-routine bratuinit.m, the definitions of PDE coefficients and Jacobian (see §3.1.3), 
and some selected commands (see bratucmds .m) to run pde2path. See also the files for detailed comments. 

(a) (b) (c) (d) 




-0-5 -0.5 -0.5 -0.5 -0.5 -0.5 



Figure 2: (a) Elementary bifurcation diagram for (18) (||ix||i;,2 over A) generated by pde2path, over a uniform 
mesh with 800 triangles. Thick lines indicate stable (parts of) branches, thin lines unstable branches, o 
bifurcation points, (b), (d) Some solution plots, (c) Preview of mesh refinement. See §3.1.5 for the quality 
of the mesh at the "ends" of branches r and the due mesh-refinement. 

3.1.3 The PDE— coefficients and Jacobians 

The coefficient function p . f is fundamental, and the jacobian function p . j ac is recommended. The 
input argument u of both is the vector of nodal values, withi/,i(-) =u(l:p.np), U2{-) =u(p .np+1 : 2*p .np) , 
. . . , un{-) =u((p.neq-l)*p.np+l :p.neq*p.np) . For the outputs c, a, /, h (of p.f ) we allow two 
forms, i.e., (arrays of) constants, or (arrays of) values on the triangle midpoints of the FEM-mesh, 
essentially as explained in [32]. There are two major ways to generate c^a^f^b from u. We first 
focus on /: 

a) Use u=pdeintrp(p.points,p.tria,u) to first interpolate u to the triangle values (again 
called u), which yields a matrix 

(Uii Ui2 ... Ui^nt\ 
UNI Un2 • • • UNnt / 

where rit is the number of triangles in the mesh. Then write, e.g., f{u) in a standard Mat lab 
way, i.e., f =-10* (u-lam*exp(u) ) ; see bratuf in Table 5. 
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b) First express, e.g., c, / as "Matlab text expression in x,y,u,ux,uy" (from [32]), with obvious 
meaning. For this, a parameter lam must be converted to a string. Afterwards, pdetxpd is 
called to evaluate the text expression on the triangle midpoints. See file bratuft.m for an 
example. 

Option a) has the advantage that it is more "natural", more flexible, and, at least for sim- 
ple expressions, slightly shorter. If, however, / depends on a::,^^;,..., then option b) might be 
shorter. To some extent it is a matter of taste which way to generate / is preferred (and similarly 
c,. . . ), therefore for (18) we provide both as templates. However, b) only allows local dependence 
f{u, ' ' ') = f{u{x)^ . . .), while in a) we have more flexibility and, for instance, can also call external 
functions in a simple way. 

For c, which in principle is the A/'xA/'x2x2 = 1x1x2x2 tensor cuki = we simply 

write c=l which corresponds to the simplest symmetric case. 

The tensor h = hij^^ i^j = 1, . . . , A/", = 1, 2 in (1) is not part of the pdetoolbox. Its coding 
and storage mimics that of c. In detail, b is an 2N'^ x m array, where m = 1 (constant case) or 
m =p.nt, in the order 

b = [hn] &112; &211; &212; . . . ^ATii; bm2] &121; ^122; • • • ; &Ar2i; bN22] • • • bim] • • • ; bNN2] , (20) 

i.e., bijk is in row 2N{j — 1) + 2i + — 2. Unlike the pdetoolbox setup for c, a there are currently no 
schemes to encode special situations such as, e.g., b (g) \/u = adxU which corresponds to advection 
into direction x. Thus, in this case, for N = 2, h reads b = [a;0;0;0; 0;0;a;0], while, e.g., ^dyU 
yields b = [0; /3; 0; 0; 0; 0; 0; Again we remark that advective terms can also be put into / such 
that p.f may simply return b=0. However, b is needed if there are advective terms and we want 
explicit Jacobians as explained next. 

If p. jsw< 3 then the user must also set p. jac to a function handle with outputs c, f u, f lam, b 
again defined on triangle midpoints, where in fact c is the same as in p.f . Compare Table 5, and 
Remark 3.1 - also concerning the meaning of fu and b. Providing the same c (and also often 
the same b) twice (in p.f and p. jac) is redundant, and there are situations where only either 
c,fu,b or flam are needed, namely jsw=l resp. jsw=2. However, we found the small overhead of 
recomputing c,b, resp. unnessessarily computing fu,b resp. flam acceptable to have a clear code. 
On the other hand, splitting the calculation of PDE coefficients and Jacobian coefficients into two 
routines is reasonable since often at least for testing it is convenient to use j sw=3 where analytical 
jacobians are never needed. 

Remark 3.1. Given coefficients c, a, /, & = c{uo)^a{uo)^ f{uo)^b{uo)^ the FEM transforms (1) into 

the algebraic system K(uo)u — F(uo) = where K is called the stiffness matrix, assembled from c, a 
and and F is the FEM representation of /. Thus, u solves the FEM discretization of (1) if the 

residual r = resi(p,u, lam) := K(u)u — F(u) = 0. The basic pdetoolbox routine to assemble K=Ko 
(in case 6 = 0) and F is [K,F]=assempde( . . . ,c,a,f ), where . . " stands for boundary conditions 
and mesh-data. To assemble the advection matrix B we additionally provide B=assemadv(p, t ,b) . 
The full system-matrix then is K = Kq — B . 

^^See, e.g., §3.5. 

^^There are various special coding schemes for diagonal or symmetric cases, see [32]. For convenience here we just 
note that in the general case c is a 4A/'^ x p.nt matrix (or in case of constant coefficients a 4A/'^ column vector) with 
djki in row 4A^(j — 1) + 4i + 2/ + /c — 6, i^j = 1, . . . , A^, k^l = 1, 2. In this scheme we would that have c=[l;0;0;l] 
to encode the Laplacian. Similarly, there are special storage schemes for symmetric a, but in general a = aij, 
i, J = 1, . . . , TV, is stored as a x p.nt matrix (resp. a A/"^ column vector in case of constant coefficients) with dij 
in row N{j — 1) -\- i. Note the somewhat non-lexicographical order. 
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Thus, if a = 6 = and c does not depend on {u, A), then with the local derivatives fu = fu and 
flam = fx returned from p. jac, the Jacobian Gu and the derivative Gx can be obtained from 

[Gu, Glam] = assempde(. . . , c, — f u, —flam), (21) 

and this is done for jsw=0. If a 7^ is independent of u (and still & = 0), then p. jac must return 

f u = - a. 

If 6 = but, e.g., a depends on then the formulas must be adapted accordingly. In other words, 
in calculating fu assume that a = in (1), i.e., if a inp.f, then define f(u) = f{u) — an and 
set fu = fu in p. jac. If c = c{u) 01 h — h{u) ^ 0, then Gu can still be assembled using c{u) and 
suitable b and fu, see §3.4 and §4.1 for examples. Similarly, -flam must always be understood in 
a "generalized" sense, i.e., it must assemble to Ga, even if this involves high derivatives of u which 
originally were implemented via c. In any case, remember that the notations fu, flam and b in 
p. jac are only conventions^ motivated by the fact that the case that only / depends on (?i. A) is 
the most common. 

For jsw=l, Gu is still assembled but Glam is calculated by finite differences.^^ For jsw = 2 
we use numjac to calculate Gu; to be efficient this requires a sparsity structure S, and here we 
assume that depends only on (all components of) u on the i-th node and all neighboring nodes, 
which corresponds to the sparsity structure obtained by [Gu,Glam] =assempde( . . . ,0,a,0) with 
aij — 1 for z, j = 1, . . . , A/^. Our experience is that numerical Jacobians are fast enough for moderate 
size problems, i.e., for up to a few thousand degrees of freedom. Of course this also depends on 
the structure of the problem: diagonal diffusion or not, weak or strong coupling of the different 
components of u. Still, assembling Jacobians is usually much faster. For jsw = 3, both Gu and 
Glam are approximated by finite differences. In any case, for both (jsw < 1 and jsw > 2) we 
assume local dependence of / on u. See, however, §3.5 for some modifications for the case of global 
coupling. J 

Remark 3.2. The boundary conditions, see §3.1.4, are updated from bc=p.bcf (p,u,lam) before 
assembling. In the (frequent) case that the BC do not change during continuation we set may 
p.bcf =@(p,u,lam) be in the init-routine (after generating be). See, however, §3.3 for examples 
with A-dependent BC. J 

As mentioned, when applicable, assembling Gu via c^a^fu and h (p.jsw=0,l) gives a matrix 
Gu^a and is faster (by orders of magnitude for large Nup) than numerical differentiation by numjac 
(p . j sw=2 , 3) , which gives a matrix Gu^n which is in general close to but not equal to Gu^a- Intuitively 
we might also expect Gu^a to be "more accurate" than Gu^n- However, we need some caution: in 
fact, Gu^n often gives better convergence of the Newton loop for the algebraic system r{u) = 

K{u)u — F{u) = 0. The reason is that Gu^a involves interpolation of nodal values to triangle 
values in c, a,6 and fu-, while for Gu^n this is done on c^a^f^b^ consistent with the definition of 
r{u). This effect becomes prominent on poor (underresolved) meshes, where the relative error 
e = \\Gu,n — Gu,a\\/\\Gu,n\\ cau be of order 0.05 or larger. However, e -> for mesh spacing 
h ^ 0. For convenience we provide the function [Gua,Gun] =jaccheck(p) which returns Gu,a^ Gu,m 
produces spy-plots of these matrices, and prints the timing and some diagnostics. 

Thus, for small Up it might appear that Gu,n is favorable. On the other hand, Gu,a is obtained 
in acceptable time on much finer meshes, where the FEM solution u^, should be much closer to a 
PDE solution u. In fact, using jsw=2,3 can even be dangerous in the sense that it may mask the 
fact that a FEM solution is not close to a PDE solution u. See §3.1.5. 

Since approximating Gx by finite differences only takes one additional call of resi, speed is not an issue for 
choosing between jsw=0 and jsw=l. In the latter case, simply set flain=0 in p. jac. See §3.2.1 for an example. 
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3.1.4 The geometry and the boundary conditions 

The domain ft is typically described as a polygon. As the pdetoolbox syntax is somewhat unhandy 
and the rectangular case is quite common we provide the function geo=rec (Ix, ly) which yields 
^ = [—Ix^lx] X [~ly^ly]' An extension is the function polygong, see [24] for its syntax. Setting up 
"arbitrary complicated" geometries is most convenient if there is a drawing img. jpg of ft in the 
current directory. Type im=imr ead ( ^ img . j pg O ; f igur e ( 1 ) ; image ( im) ; [x , y] =ginput ; which 
yields a crosshair. Click (counterclockwise) on dft^ stop with return. The obtained vectors x,y 
can be saved as a *.mat or *.txt file, and piped through geo=polygong(x,y) . Finally, geometries 
can also be exported from the pdetoolbox GUI. 

The pdetoolbox syntax for the boundary conditions (2) is also somewhat unhandy. For 
scalar equations the most common BC are homogeneous Dirichlet or Neumann BC. The routine 
[geo,bc]=recdbcl(lx,ly,qs) approximates Dirichlet BC over rectangles via Robin BC of the 
form n • (c Vu) + qghu = with a large = qs. 

For c of order 1 typically Qs = O(IO^) or Qs = O(IO^) works well. For homogeneous Neumann 
BC we provide [geo,bc]=recnbcl(lx,ly), with the extension [geo,bc]=recnbc2(lx,ly) to 2- 
component systems. 

For the genuine systems case (or the case of non-rectangular domains) we provide the routines 
bc=gnbc(pneq, varargin) and bc=gnbcs(pneq, varargin) . For a system with neq components 
and a domain with nedges edges, bc=gnbc(neq,nedges,q,g) creates "generalized Neumann BC" 
(2) that are given as numerical data. Different boundary conditions qj^gj at the edge with index 
j are generated by a call of the form bc=gnbc (neq 

> 9.1 'Si' • • • 5 9.nedges 3 gnedges ) • FoT q givCU in 

terms of an explicit formula, e.g., involving x^u^Vu^ the function gnbcs accepts a string variable 
encoding of g and q or gj, q^- and otherwise works in the same way. 

3.1.5 Error estimates and mesh adaption 

As an ad hoc way to check whether a FEM solution =p.u approximates a PDE solution u 
we provide [q,ud] =ineshcheck(p, cmp) . This (adaptively) refines the FEM mesh in p to roughly 
the double number of triangles and calculates a new FEM solution Uh^new from the old solution 
'^/i,oid- Then Uh^oid is interpolated to Uh,new on the new mesh, iXdiff — '^/i,new — Uh^new is formed, 
and ||i/diff llcx) and the relative error ||t^diff ||oo/||'^/i,new||cx) are printed. The new solution structure q 
and the difference ud= u^i^ are returned, and for cmp> we additionally generate a plot of the 
cmp*^ component of t^diff- For instance, in Fig. 3a) we check the mesh of point 20 on the q branch 
in Fig. 2 which indicates that the mesh in 2 corners is clearly too poor, and before continuing for 
smaller A we should refine the mesh. 

Thus, some (automatic) mesh adaption may be vital for reliable continuation. The pdetoolbox 
comes with mesh refinement based on an a posteriori error estimator as follows. For the scalar 
Poisson problem —Au = /, u\oq = 0, let be the FEM solution and u the PDE solution. Then, 
with Q^,/3 > some constants independent of the mesh, 

||V(« - Uh)\\L2 < a\\hf\\ + PDhiuh), (22) 

-| /<-) 

where h = h{x) is the local mesh size, and Dh{v) — (X^^^^. h'^{dn^v)'^) . Here dn^v is the jump 
in normal derivative of v over the edge r, hr the length of the edge, and Ei the set of all interior 

^^See bratucmds.m for the calling sequence for Fig. 3a). It is also often useful to repeat calls to meshcheck, 
i.e., here call next q=meshcheck(q, 1) , until the relative error becomes suffiently small, indicating that Uh is a good 
approximation. 
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edges. For equations — V(c Vu) + au — f this suggests the error indicator function 



1/2 



E{K) = a\\h{f - au)\\K + P 



redK 



cVuhf 



(23) 



for each triangle, which is calculated by the Matlab routine pdejmps. For convenience we provide 
the interface routine err=errcheck(p) .^^ Calling, e.g., err=errcheck(p) yields err=0.273 which 
somewhat overestimates the error in Fig. 3(a). In cont, for errchecksw>0 we call errcheck after 
each successful step and store the result in p. err. 

The (basic) mesh refinement strategy then is to introduce new triangles where E{K) is large. 
This is done by the pdetoolbox routine ref inemesh, but we provide the interface routine p=meshref (p, 
varargin).^^ Since meshref also interpolates the tangent r to the new mesh, we can continue 
with cont immediately after mesh refinement. However, instead of mesh refinement, which means 
introduction of new points into the mesh, we rather need mesh adaption, which means refine- 
ment where necessary, but coarsening where possible, to limit storage requirements. In pde2path, 
mesh-adaption is implemented in an ad hoc way in the function p=meshadac(p, varargin) by first 
interpolating a given solution to a (typically somewhat coarse) "base-mesh" or "background-mesh" 
and then refining. 

During continuation runs there are basically two strategies for this, which can also be mixed: 

(i) call meshadac every p.amod*^ step, for p.amod>0, or 

(ii) call meshadac whenever p . err>p . errbound (choose p. errchecksw>l for this)-^^. 
See Figures 3(b), (c) for a comparison of the different approaches. 

(a) 

meshcheck(qO,1) 





Figure 3: (a) Error in u at qO=point 20 on the q branch from Fig. 2 obtained from calling meshcheckCqO , 1) . 
(b) p. err over A for continuing from point 10 on the q branch, without mesh- adapt ion (labeled qO), with 
strategy (i) (amod=10, labeled ql) and with strategy (ii) (p . errbound=0 . 1, labeled q2). (c) The bifurcation 
diagrams belonging to (b). In ql a rather large jump appears from refinement after the 10^^ step. 

In summary, mesh adaption strategies and error bounds are highly problem dependent, and 
moreover, may not be rigorously justified for the system case or general BC. Thus, although it 

With a,f ,b returned from p.f, this also re-calculates / (via f b=bgradu2f (p,f ,b,u) ) to include 6 Vi^ into 
fb, as pdejmps does not take 6 Vix into account directly. Similarly, the mesh refinement below always recalculates 
/ to include h Vu. Moreover errcheck also contains our settings for some tunable parameters of pdejmps. 

^^It is not a priori clear if this is also suitable for systems. Nevertheless we found pdejmps to work well also for 
the problems considered below. 

^^where varargin takes pairs ^maxt^ ,maxt=number of triangles aimed at, or ^ngen^ ,ngen=number of refinement 
steps, or ^ eb ' , eb=error bound. Calling for instance qOr=meshref (qO, ^eb^ ,0.0025, 'maxt' , 50000, ^ngen^,20) 
shows that to achieve an estimated error< 0.0025 we need about 30.000 triangles. 

^^base-mesh given by p.bpoints, p. hedges, p.btria, varargin as above. The implementation of some true 
adaption routine is on our to-do- list. 

where we refine until err<p . errbound/2 in order to allow some margin for the next steps 
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works well for a number of examples we considered, mesh-adaption should be applied with care. 
Also note that mesh-adaption may lead to spurious jumps in detA. 



Remark 3.3. For bifurcation from trivial branches, another good strategy is to prepare a/inerbase 
mesh than the starting mesh, for instance if the trivial branch consists of spatially homogeneous 
solutions, but the bifurcating solutions develop sharp gradients. For convenience we also provide 
the functions p=newinesh(p) , which interpolates the current (u^r) to a new mesh generated after 
user input in the form hmax or nx,ny, and p=setbmesh(p) which sets the base mesh to the current 
mesh. This should be called if it is expected that the current mesh is a good base for adaption in 
the steps to come. J 

3.1.6 The linear system solvers 

Recall that after discretization with rip points we have nodal values u E MP with p = Nrip large, 
and 



are large, but sparse (block) matrices. The question is how to best solve GuV = t and the bordered 
systems such as At — z, respectively. In all the examples that we considered, our experience is 
that the highly optimized matlab solver z=A\b of Az — b works remarkably well, but for easy 
customization of the code we never call \ directly but use two interface routines: 

1. v=p.lss(M,r,p,lam) to solve Mv = r with M = G 

2. z=p.blss(A,b,p,lain) to solve Az = b with A e W^^^'p^^. 

Here blss and Iss stand for (bordered) linear system solver. 

The default solvers Iss and blss just contain one command, namely v=M\r resp. z=A\b. Nev- 
ertheless, for large systems or for some special classes of problems iterative solvers might work 
better, and as templates we provide the two routines ilss and iblss, using gmres with (incom- 
plete) LU factorization luinc as preconditioners. These should, of course, be reused as long as 
gmres converges quickly, and here (and in resinj.m) we thus introduce some global variables, 
namely global L U; resp. global bL bU. Thus, when using, e.g., ilss the user must also issue 
global L U; L=[] ; U=[] ; from the command line. 

It turns out that for scalar problems these outperform the direct solvers Iss and blss for large 
> 10^, say. On the other hand, for systems, luinc becomes exceedingly slow such that 
\ beats gmres with LU preconditioning even for very large Up. In summary, the iterative solvers 
ilss and iblss should only be regarded as template files to create problem specific iterative solvers 
when needed. See also §3.5 for adaptions of Iss and blss to some special situation. 

Finally, various approaches have been proposed for the solutions of the bordered systems Az = 6, 
see, e.g., [16, 12]. As alternatives to blss we provide bellss ("bordered elimination") and belpolss 
("bordered elimination plus one"). To use these, simply set, e.g., p.blss=@belpolss. In our tests 
the performance of bellss and belpolss is roughly the same as blss. 

3.1.7 Screen output, plotting, convergence failure, auxiliary functions 

The screen output during runs is controlled by the two functions p.headfu (headline) and the 
function p.ufu. These are preset in stanparam as p.headfu=@stanheadfu, p.ufu=(§stanufu to 

^°The reason for this construction is that we do not want to make L, [/ a part of p since this needs a lot of disk 
space when saving p: typically, we get fill-in factors for L,U of 10 and larger. 



Gu e MP^f and 




(24) 
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first print a headline and then, after each step, some useful information. To print some other 
information the user should adapt stanheadfu and stanufu to a local copy, say myhead.m, and set 
p.headfu=(§myhead, and similar for stanufu and p.ufu.^^ The bifurcation diagram and solution 
plots are also generated during continuation runs, but in general it is more convenient to postprocess 
via plotbra, plotsolf etc. 

The files p* .mat and bp* .mat contain the complete data of the respective point on a branch. 
Thus, a run which is no longer in memory can be simply reloaded by, e.g., q=loadp(pre,pname, 'q'), 
where pre , pname is the name data of a previously saved point, and the third argument is used to 
set the directory name for the newly created struct. The loaded point will often be either the last one 
or the first; in the latter case, to change direction of the branch, use, e.g., q=loadp( 'p' , 'pi ' , 'qO ; 
q.ds=-q.ds; q=cont(q); 

If the Newton-loop does not converge even after reducing ds to dsmin then cfail.m is called. 
The standard option is to simply abort cont, but we offer a number of alternatives, e.g., to change 
some parameters like dsmin or imax, or to try, e.g., some mesh refinement or adaption. Clearly, the 
choice here is strongly problem dependent, and thus we recommend to adapt cfail.m if needed; 
see §3.6 for remarks on such "customization without function handles". 

Besides those already mentioned we provide further auxiliary functions, see [23, m2html] for a 
complete documented list. 

3.2 The Allen— Cahn equation with Dirichlet boundary conditions (ac) 

In our second example we use Dirichlet boundary conditions (DBC), and explain some ad hoc 
parameter switching, and time integration. We consider a cubic-quintic (to have folds) Allen- 
Cahn equation 

—jj.Au — Xu — ^ = on = [—Lx, L^] x [— Ly, Ly], u\q^ 0, (25) 

with two parameters fi > and A G M. We use [p .geo ,bc] =recdbcl (Ix, ly , le3) to approximate 
the DBC, and set L^^ = 1 and Ly = 0.9 to break the square symmetry present in bratu in order to 
have only simple bifurcations, namely at X^i = iJiTi^{{k/ + {l/LyY). First we fix /i = 0.25 which 
yields An = 1.3784, A21 = 3.2289, A12 = 3.6630, . . ., and continue in A, which yields Fig. 4(a)-(c). 
After branch switching we turn on mesh-adaption after each 5 steps. See acdemo.m or accmds.m 
for more details, which also contain an example of perturbing a solution and subsequent time 
integration. 



(a) (b) (c) (d) 




Figure 4: (a) Elementary bifurcation diagram for (25) with /i = 0.25. No secondary bifurcations occur, and 
the mode-structure on each branch is completely determined at bifurcation. (b),(c) some points on branches 
as indicated, (d) Solution after continuation in ji from (b) to /i = 0.1. 



^^For p.timesw>0 we also plot some timing information at the end of cont. 
Obviously this is often quite redundant, but it is necessary if, e.g., some mesh refinement occured during 
continuation. To save disk space, however, we deliberately chose to not make Gu a part of p. See also footnote 20. 
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3.2.1 Parameter switching 

Unlike AUTO, pde2path (currently) has no switches or presets for multi-parameter continuation. 
However, switching to a new parameter for continuation can be achieved in a simple and flexible 
way by modifying the structure p from the command line. As an example we want to continue in /i 
from point 10 on q to /i = 0.05. This is achieved by the commands in Table 6, and yields Fig. 4(d). 
The basic idea is to copy q to w (this is not strictly necessary) and then reset w.f and w. jac. We 
make our life simple by setting w. jsw=l such that we do not need and set ^ = 10~^ since the 
dependence on /i is quite sensitive. 

We also introduce a new parameter w.upl ( 'user parameter 1', but any name will be fine) which 
is used to pass the current A to acf mu, see Table 6. 

w=q;w.upl=w.lain;w.lam=0.25;w.laiiimin=0.05;w=setfn(w) ; w.ds=-0.01; 

w . f =@acf mu ; w . f =@ac j acmu ; w . j sw=l ; w . parasw=0 ; w . xi=le-6 ; w . restart=l ; w=cont (w) ; 

function [c ,a,f ,b] =acfmu(p,u,lam) % AC for cont in mu 
u=pdeintrp(p.points,p.tria,u) ; c=lam; a=0; f =p.upl*u+u. "3-u. "5; b=0; 

function [c,fu,f lam, b]=acj acmu (p,u, lam) % Jacobian for cont of AC in mu; run with jsw=l 
u=pdeintrp(p.points,p.tria,u) ; c=lam; flam=0; b=0; f u=p.upl+3*u. "3-5*u. "4; 

Table 6: Switching to continuation in /i, commands, and modified coefficient and Jacobian functions. 

3.2.2 Time integration 

For time integration of (3) using the struct p we provide a simple semi-implicit Euler method. 
Writing u^'^^ for u{tn^ choosing a time-step /i, approximating dtu{tn) ^ j^{u^'^^^^ — u^'^^) where 
= ^0 + and evaluating, e.g., V • (c (8) Vu) as V • [c{u^'^^) (g) Vix(^+^)) we obtain, on the FEM 
level, 

iM(^(^+^)-^(^)) = -i^(^(^))^(^+i)+F(^(^)) 
^ = (M+/ii^(n(^)))-i(Mt^(^)+/iF(n(^))). 

Here M is the mass matrix and K is the stiffness matrix on time-slice n. This is implemented 
in tint (p,h,nstep,pmod) , where nstep is the number of time steps, and a plot (of component 
p.pcmp) is generated each pmod'th step. Thus we may call, e.g., p.u=p.u+0 . l*rand(p.neq*p.np, 1) ; 
p=tint(p,0. 1,50,4,10) to first perturb a given solution and then time-step. See accmds.m for 
an example, where we perturb a solution on the unstable part of the q branch into both directions 
of the unstable manifold; in the subsequent time integration the solution converges to the stable 
trivial solution or the stable q branch, respectively, as expected. However, the main purpose of 
tint is to generate (stable) initial data for continuation, i.e., after tint call cont. See also §5.2 
for an example where tint is used in this spirit. 

Remark 3.4. In tint we assemble ^(m^")) and solve (M + = y^") in each step 

by Iss. Clearly, for special cases this can be optimized: for instance, if c, a,6 do not depend on 
then the textbook approach would be to assemble K at the start, followed by some incomplete 
LU-decomposition of M + hK combined with some iterative solver. However, similar remarks as 
in §3.1.6 apply, and thus we use the very elementary form above, but stress again that tint in its 
present form is not intended for heavy time-integration. J 
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3.3 The Allen— Cahn equation with mixed A-dependent boundary conditions 
(achex) 

We illustrate a few more possibilities with pde2path by modifying the Allen-Cahn example (25) 
from above. We consider again (25), i.e., —0.25Au — Xu — -\- = 0, but instead of homogeneous 
Dirichlet BC on a rectangle, we consider hexagonal domains and parameter dependent mixed 
Dirichlet-Neumann BC. Figure 5(a) shows an example for fi, which basically consists of a square, 
with the top boundary shifted by Sy = 0.5 between [—IxJx]^ Ix — 0.5. Denote this part of by 
Td, and set Tn = dQ\rD. To define the domain we could, e.g., use the pdetool GUI to draw a 

(a) (b) (c) (d) 




Figure 5: (a) A hexagonal domain, already with a basic mesh, created by p=staninesh(p,0 . 1). (b) Bifur- 
cation diagram for (25) with BC given by (26). (c),(d) selected solution plots. 

polygon composed of six edges, one for each segment and export the geometry. However, usually 
the function polygong [24] is much more convenient. Second, we want to define the boundary 
conditions 

n • Vu = on Fat, u = Xx on Tjj. (26) 
To implement this we use a stiff spring approximation on F^) in via gnbcs, i.e., 

qd=inat2str ( 10^4) ; gd= [inat2str ( 10^4*lam) ' *x ' ] ; qn= ' ' ; gn= ' ' ; 
bc=gnbcs(l,qn,gn,qn,gn,qn,gn,qn,gn,qd,gd,qn,gn) ; ^ ^ 

With pde2path we perform a continuation starting from the trivial zero solution and obtain the 
bifurcation diagram plotted in Fig. 5; see ac6cmds. Bifurcation detection and branch switch- 
ing work without problems, and the error estimate is always well below 0.01. To generate both 
parts of the r branch we first call rl=swibraCp' , 'bp2' , 'rl' ,-0. 1) ;rl=cont (rl) and then 
r2=loadpCrl ' , 'pi ' , 'r2' ) ; r2.ds=-rl .ds; r2=cont(r2) to proceed in the other direction. At 
the end of acGcmds we also run an example with = A on F^ implemented via gnbc. 

3.4 A quasilinear Allen— Cahn equation (acql) 

To give an example of a more complicated Jacobian we modify (25) to the quasilinear Allen-Cahn 
equation 

-V '[{0.25 + 5u + -fu^)Vu]- f{u,X) = on Q = [-L^, L^] x [-Ly, Ly], u\dn = 0, (28) 

with f{u, A) = Xu + u^ — u^ and Lx — l^Ly — 0.9 as before. See acqlf .m. The linearization around 
u gives the linear operator 

Gu{u,X)v = -V-[(0.25 + 5ix+7ix2)Vi;] + [-/^(ix, A)-5Aix-27(Vi/-Vi/+i/Ai/)]i;-[(5+27i/)Vix]-Vi;. 

^^see geo=hexgeo(lx,dely) , which also contains a shghtly edited output of the GUI for comparison. 
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Hence, in acqljac .m we now have f u = /^^ + 5Au + 2'^{Vu • Vu + uAu), and bm = (5 + 2^u)ux 
and bii2 = (5 + 2^u)uy, cf. Remark 3.1. To generate (ux^Uy) and Ai/, as coefficients in acqljac. m 
we use pdegrad resp. pdegrad, pdeprtni and pdegrad again, see [32]. 

The term 6u in c changes the u i-^ —u symmetry of the Ahen-Cahn equation (25). The 
bifurcation points from the trivial branch u = in (28) are as in (25), but the bifurcations change, 
see Fig. 6. In particular the first bifurcation changes from pitchfork to transcritical. 



u at wpIO 




Figure 6: Elementary bifurcation diagram for (28) with S = —0.2 and 7 = 0.05, and some solution plots. 
The two blue branches are in fact one branch, and the corner at the transcritical bifurcation from the trivial 
branch is due to the choice of vertical axis. The symmetry u —u no longer holds, and "up" humps are 
steeper than "down" humps due to ^ < 0. The w-branch is still double due to the x —x symmetry. 

3.5 An Allen— Cahn equation with global coupling (acgc) 

As an example of a "non-standard" elliptic equation we treat an Allen-Cahn equation with a 
global coupling. We fix /i = 0.1 and A = 1 in (25), introduce a new parameter (again called A) and 
consider 

G{u, A) := -O.lAu - u - + - X {u) = on fl^ [-vr/2, v^/2]^ uIqq = 0, (29) 

where {u) = J^udx. The term X{u) is called a global coupling or global feedback, positive for 
A > resp. negative for A < 0. Problems with global coupling occur, e.g., in surface catalysis, where 
global coupling arises through the gas phase [25], in semi-conductors and gas-discharges [34, 31], 
and as "shadow systems" in pattern formation when there is a very fast inhibutor diffusion [15]. 

The global feedback does not fit into the framework of (1) if / is assumed to be local. For the 
definition of G(u) this is not yet a problem as we may simply define / as, e.g, 

f =u+u. ''S-u. ''S+lamH^triint (u,p.points,p.tria) ; 

where triintCg, points, tria) is the Riemann sum of J g{x)dx over the given mesh. However, 
for continuation we make extensive use of Jacobians, and Gu(u) is now given by 

[Gu{u)v]{x) = -0.1Av{x) - (1 + 3u{xf - 5u'^{x))v{x) -X{v). 

As yet we cannot deal with last term, cf. Remark 3.1. The first try would be to simply ignore it in 
continuation, but this in general only works for small |A| while for larger |A| we loose convergence 
in the (false) Newton loop. We can express {v) on the FEM level via a matrix M such that 
Gu{u)v ^ [K — XM)v. Essentially, for "natural parametrization" we need to solve 

Gu{u)v = r, where Gu{u) ^ {K - Xvr]^) with z/, 7/ G . (30) 
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Here 77 = {oT)^ where (ai, . . . , J contains the triangle areas, T G W^^^'^p interpolates u G W^^ 
from nodal values to triangle values (such that {u) =triint (g, points ,tria)=eta*u), and vi — 
l(f)i dx corresponds to adding (v) to all nodes with the correct weight. However, {K — Xur]^) is a 
full matrix and should never even be formed. Instead we customize Iss to use a Sherman-Morrison 
formula which gives (for (30)) 

V = K-^r + a(K-^u)(r]^K-y, a = ^ . (31) 

1 — AT]^ K 

In acgcjac.m we then just ignore the term A (ix). Similar remarks apply to the bordered systems 
solved by blss. 

In the actual implementation we introduce global variables nu,eta. The idea is that it is 
sufficient to calculate once for a given mesh, for instance in acgcf .m, as this is always called 
before the Jacobian acgcjac or the linear system solvers gclss or gcblss. If we set aside mesh- 
refinement then we could calculate z/, 77 at startup and store them e.g. as p . nu , p . eta but with mesh 
refinement global variables are more convenient. See Table 7 for the full code for acgcf .m, Iss.m 
for this example, and Fig. 7 for the result of the basic continuation runs contained in acgccmds .m. 
We switch off spectral calculations and bifurcation checks by setting spcalcsw=0; bif checksw=0; 
since out of the box these would be based on the (wrong) local Jacobian, and the two branches 
were generated by using two different starting points. 

function [c,a,f ,b] =acgcf (p,u,lam) % AC global coupling 
global eta nu; try se=size(eta,2) ; catch; eta=[]; se=0; end 
if (se'^=size (u, 1) ) % eta not yet set, or mesh is refined 

C=n2triainat (p.points,p.tria) ; ta=triar(p. points, p. tria) ; eta=ta*C; 

[M , nu] =asseinpde (p . be , p . points , p . edges , p. tria, 0,0,1) ; end 
uin=eta*u; u=pdeintrp(p.points,p.tria,u) ; c=0.1; a=0; b=0; f =u+u. "3-u. "5+lain*uin; 

function x=gclss(A,b,p,lam) % Iss for AC with global coupling, Sherman-Morrison 
global eta nu; y=A\b; z=A\nu; alO=lam*eta*z ; al=lam*eta*y/(l-alO) ; x=y+al*z; 

Table 7: Definition of / and customized gclss. m; see also acgcjac.m and gcblss. m. 



(a) (b) (c) (d) (e) 




Figure 7: (a) Two solution branches for (29), and four selected solutions. By positive global feedback, the 
plateau in (b) {u around 1.93) is substancially above the zero (1 + V^)/2 ^ 1.62 of f{u) = u -\- — u^. 
Here, some mesh refinement near the boundary is also crucial. Decreasing A to slightly negative values u 
gets pushed below near the boundary (c). On the other branch some somewhat localized solutions are 
found (d), (e). Note that (29) is symmetric w.r.t. {u^X) (— li. A). 
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3.6 First summary, and some remarks on customization 

We end this introductory section based on scalar examples with a first summary and some imple- 
mentation remarks. 

The p . f =0 . . syntax has the advantage that multiple version of f can be maintained and switch- 
ing can be done by only changing p . f =@ . . . On the other hand, we do not want to overwhelm the 
user with such options, and thus we restricted the "user-definable" functions to p . f , . . . ,p .headf u 
from Table 2, where in fact in most cases the user only needs to set up p . f , and p . jac for p . j sw< 2. 
Nevertheless, as outlined above any function of pde2path can be customized for a given problem 
by just copying it from . . /p2plib/ to the current directory (where Mat lab searches first) and then 
modifying it. Main candidates for customization are, e.g., plotbra.m, plotsol.m if additional 
features/options are desired in plotting the bifurcation diagram or/and the solutions. See, e.g., 
§5.1 and §5.2. 

Most functions of pde2path only require a few input/output arguments. An important excep- 
tion is plotbra(p,wnr,cmp,varargin) where varargin is a possibly long list of argument /value 
pairs. See plotbra.m for a detailed description, and also the various plotbra example calls in the 
demos. 

As mentioned, by default there are no global variables in pde2path, with the exceptions p j , lamj 
which are set for numerical differentiation in resinj, and possibly LU preconditioners for iterative 
linear system solvers, see §3.1.6. On the other hand, e.g., p.f, p. jac, p.lss etc. do not return 
p. This is to have a somewhat clean distinction between functions for specific calculations and 
function like p=cont (p) , p=swibra( . . . ) , p=meshref (p) which modify the structure p, including 
the mesh. As a result, the user might want to introduce some global variables to streamline 
calculations, see §3.5 for an example. These should then be declared before initialization of p. 

For convenience, in Table 8 we summarize the typical steps in the usage of the software. 

Initialization. Declare (user defined) global variables (if any). Initialize structure p, typically by first 
calling p=stanparain(p) , followed by problem dependent calls to define (function handles for the) PDE 
coefficients, BC and Jacobian, and the geometry, mesh, and starting point, 
function p=cont(p) 

1. If restart=l then inistep: generate first two points on branch and (secant) tq 

2. Predictor (u^^X^) = {uq^ Aq) + dsro with stepsize ds 

3. Corrector: depending on parasw and Aq use nlooppde for (16) or nloopext for (10) or (13). 
This uses getder , getGu resp. getGlam to obtain derivatives, and Iss resp. blss as linear systems 
solver. 

4. Call sscontrol to assess convergence (res, iter returned from nlooppde resp. nloopext): If 
res<p.tol accept step, i.e., goto 5, (and increase ds if iter<imax/2). If res>p.tol and ds> 
dsmin then decrease ds and goto 2. If res>p.tol and ds=dsmin then no convergence, hence call 
cf ail. 

5. Postprocessing: calculate new tangent ri by (12), call spcalc (if spcalcsw=l), bifdetec (if 
bif checksw=l). Check for error and mesh adaption. Update p, i.e., put uq = ui, Xq = Xi^ 
To = n into p, call out=outfu(p,u,lam) , plot and save to disk. Call p.ufu for printout and 
further user-defined actions. 

6. If stopping criteria met (p.ufu returned 1 or stepcounter>nsteps) then stop, else next step, i.e., 
goto 2. 



Post— processing. Plot bifurcation diagrams via plotbra (plotbraf) and solutions via plot sol 
(plotsolf ). If bifurcations have been found, use swibra (and cont to follow some of these). 



Table 8: Typical software usage, including pseudo-code of p=cont(p), with main function calls. 
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4 Some prototype Reaction— Diffusion Systems 



Pattern-formation in Reaction-Diffusion Systems (RDS), in particular from mathematical biology 
[22] , is one of the main applications of path- following and bifurcation software. Here we first consider 
a quasilinear two-component system with "cross diffusion" from chemotaxis [20] to explain the setup 
of c in this rather general case, and the setup of general domains in pde2path. We essentially recover 
the bifurcation diagrams from [20] without special tricks or customization. 

Our second example is the Schnakenberg model [29], which is semilinear with a diagonal con- 
stant diffusion matrix, and thus in principle simpler than the first example. However, here we are 
interested in a more complete bifurcation picture, and the Schnakenberg model shows many bifur- 
cations already on small domains. Therefore we need some adaptions of the basic cont algorithm 
to pmcont (parallel multi continuation), and we introduce f indbif to locate some first bifurcations 
from the homogeneous branch. 



4.1 Chemotaxis 

An interesting system from chemotaxis has been analyzed in [20], including some numerical path- 
following and bifurcations using ENTWIFE. The (stationary) problem reads 

n n( \\ (D/\ui- \V • {uiVu2)\ (rui{l-ui)\ 

where A £ M is called the chemotaxis coefficient and D > and r £ R are additional parameters. 
In (32) we have 6 = 0, may set o = 0, and identify the second term with /(«). The linearization 
reads 

(33) 

and in the notation from Remark 3.1, the first matrix in (33) relates to c, the second to — fu, and 
the last term gives —b (g) \/v with bm = —XdxU2^ &112 = —XdyU2, and bijk = else. 

4.1.1 Bifurcation diagram over rectangles (chemtax) 

Following [20] we first study (32) on a rectangular domain ft = [—Lx/2, Lx/2] x [— L^/2,L^/2] 
with homogeneous Neumann BC. Again a number of results can then be obtained analytically. 
There are two trivial stationary branches, namely u = (0,0) which is always unstable, and u = 
t^* = (1,1/2). From the BC, the eigenvalue problem Mv = /iv for the linearization around 
t^* has solutions of the form /i = ii{m^l^\)^ v = v{m^l^X;x) = (^em,i{x,y) with G and 

emA^^y) = ^^^(^(^+ ^"^'^^ " (1,0), (0,1), (2,0), (1,1),.... To 
study bifurcations from we solve /i(m, / , A) = for A which yields 

A^,, := 4(L»A;2 + r){e + l)/k'', where ^ := tt^ (^ + il^ . 

As in [20, Fig.3] we choose D = 1/4, r = 1.52 and the "1 x 4" domain Lx = 1, Ly = 4, which yields 
Table 9. 

To encode (32) we note that cim = C1122 = -C*,ci2ii = C1222 = -A«i,C22ii = C2222 = 1, and all 
other entries of c are zero. In particular, c is isotropic and thus we may use isoc.m to encode it. 
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(0,2) 


(0,3) 


(0,1), 


(1,0),(0,4) 


(1,1) 


(1,2) 






12.01 


13.73 


17.55 


17.57 


18.15 


19.91 





Table 9: Bifurcation from ix* = (1, 1/2) in (32), D = 1/4, r = 1.52. 



see Table 10 for chemf .m. For convenience and illustration, here by default we first use p. jsw=3 in 
cheminit . m such that p. jac need not be set. With p=staiimesh(p, 0.075) leading to p.nt=2376 
this still gives quick results, which moreover essentially do not change under mesh refinement. 
Also, in cheminit we introduce p.vol=|fi|; we want to use this quantity in chembra.m since the 
bifurcation diagrams in [20] plot \\ui — l||^i/|fi| over A. Again this is a simple example that the 
user can augment the structure p with whatever is useful. The commands in chemcmds . m yield the 



function [c,a,f ,b] =cheinf (p,u,lain) % chemotaxis system with isoc 
u=pdeintrp(p .points ,p .tria,u) ;a=0;b=0; vl=ones (1 ,p .nt) ; 
fl=r*u(l, :) .*(l-u(l, :)) ;f2=u(l, :) ./(l+u(l, :))-u(2, :) ; f=[fl;f2] ; 
D=0.25;r=1.52; c=isoc ( [ [D*vl -lam*u(l , : )] ; [0*vl vl] ] ,p .neq,p .nt) ; 



Table 10: chemf .m as a prototype for definition of PDE coefficients in case of a (nonsymmetric) c 
depending on u and A. See isoc and the assempde documentation for the order of Cijki in c. 

bifurcation diagram in Fig. 8, where the bifurcation values Xjn,i (except for Aio = A04) are found 
with resonable accuracy, and which agrees well with [20, Fig. 3(a)], with one exception: on the (1, 1) 
branch there is a loop near A = 20.5 with two bifurcations, during which the solution structure 
changes as detailed in (b),(c). Presumably, this loop was just missed in [20] due to a larger stepsize. 




(b) 

u at q4p10 



(c) 

u at q4p25 



-0.5 0.5 



Figure 8: (a) Bifurcation diagram for (32) with jsw=3, i.e., numerical Jacobians, and nt=2376. Of the 
bifurcating branches only the (0, 2)-branch is stable in a certain A range, and a number of secondary bifur- 
cations occur on each branch. (b),(c) The shape of solutions before and after the loop on the (1, 1) branch. 
For jsw=l we need finer meshes (nt ^ 10^ to 2 • 10^ and adaptive refinement), which, while giving smaller 
error-estimates, also destroy the speed advantage of assembled Jacobians. 

Alternatively, to run (32) with jsw=l we also provide chemjac.m which encodes (33). For this, 
however, we need considerably finer meshes, mainly since the calculation of the coefficient Au2 
(needed for jsw<2) via pdegrad and pdeprtni does not go together well with Neumann bound- 
ary conditions, since the averaging involved in pdeprtni produces some error at the boundaries. 
Therefore, we also replace p=cheminit (p) by p=cheminitj (p), which resets a number of switches 
to (re)run (32) with jsw=l. See the end of chemcmds. m resp. chemdemo.m. 
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4.1.2 Drawing general domains ( animal chem) 



We now consider (32) on the animal-shaped domain in Fig. 9, taken from [22], with Neumann BC. 
To set up we proceed graphicahy as explained in §3.1.4, see animalgeo .m, also for the setup of 
the BC. The plots in Fig. 9 are generated from the commands in animal cmds .m. For problems of 
this type, the bifurcation directions from a trivial branch are often most interesting. 




(b) 



(d) 



at bp1 



at qp25 



at bp2 






1 2 3 



Figure 9: Bifurcation diagram for (32) {\\ui — over A), first bifurcation direction from the trivial 

branch, stable solution on the bifurcating branch, and the second bifurcation from trivial branch. The 
bifurcation points on the trivial branch are quite close. Thus, on the trivial branch we need to run cont 
with small dlammax. 



4.2 The Schnakenberg model (schnakenberg) 

We consider the (stationary) Schnakenberg system in the form 

= G(^)=f"i^^"^^r1^^V (34) 

We use A as bifurcation parameter, fix d = 60 and consider (34) for (x, ?/) G = [—Ix^ Ix] x [~ly^ ^y] 
with Neumann BC. Over = the spatially homogeneous solution ?i*(A) = (A, 1/A) becomes 
Turing unstable [22] when decreasing A below Ac ^ 3.2, with critical wavenumber fee ^ 0.63. In 
2D, the most famous Turing patterns are 

u{x, y) = u* ^ A cos{kcx) + h.o.t. (stripes) , 

u{x^y) = i/* + Acos{kcx) + Scos(^x) cos{^kcy) + h.o.t. (hexagons, or hexagonal spots), 

where A,B ^M? are suitable amplitudes, h.o.t. stands for higher order terms (in A,B,X — Ac), and 
we dropped the A dependence of all terms. Over = [—IxJx] x [~J^yJy]i if the domain size and 
BC permit it, both (spots and stripes) bifurcate from the trivial branch at A = Ac- Here, to make 
Ac a simple bifurcation point (for vertical stripes), we choose Ix — 2rmi/kc and ly — 2n67T / {y/3kc) , 
m^n ^ N ^ where 5 ^ 1 is a deformation parameter, such that for 5 7^ 1 the double bifurcation 
point splits into two simple bifurcation points. To locate these, we start on the homogeneous 
branch and use a bisection type routine based on the number of negative eigenvalues of Gu^ see 
f indbif .m. Figure 10 shows a bifurcation diagram and some solution plots obtained for m = n = 2 
and 5 = 0.99. 

Here the problem is that using the standard cont algorithm we quickly obtain some undesired 
branch switching. For instance when continuing the stripe branch s with standard settings we 
switch to the beans branch b when approaching its bifurcation point. This particular branch 
switching can be avoided by decreasing ^ to ^ = 0.1/p.np, say, but only to the effect that we get 
branch-switching at some later point on the s branch. Such undesired branch-switching is a serious 
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-15 15 -15 15 -15 15 -15 15 -15 15 -15 15 



Figure 10: Bifurcation diagram and a selection of patterns for (34). The branch f in the bifurcation 
diagram represents the homogeneous solutions, s the stripes, s2 the phase shift of s, h the hot spots, c the 
cold spots, b and b2 are mixed modes, also called beans, o are bifurcation points. Thick lines are stable 
and thin unstable. For s, h and b we plot the maximum of Ui and for s2, c and b2 the minimum. As 
usual, hp3 stands for the third point of h, bbp2 for the second bifurcation point of b and ri at f bp3 for the 
tangent in the third bifurcation point of f . On the right we also plot \uki\^^'^, where u is the discrete Fourier 
transform of Ui — (i^i), see four.m in directory schnakenberg. These Fourier plots are often interesting for 
pattern forming systems: for instance b2pl0 shows that the pattern is essentially still generated from the 
basic harmonics exp(i/cc^)^ and exp(i/cc(^x ± ^y)^^ m^n = ±1, while at hp24 a rather large number of 

higher harmonics exp(i/ccx)^ exp(i/cc(^x ± ^y))^, m, n G N, contribute. Branch switching at some of the 
further bifurcation points yields a number of further interesting patterns, including some "snaking" between 
stripes and hexagons, see [33]. 

problem in all continuation algorithms, see, e.g., [30, §3], and we use a modification pmcont of cont 
explained in the next section, which also incorporates some parallel computing for speedup. 

4.3 pmcont 

Theorem 4.4 in [16] guarantees that the standard continuation converges to a given branch for 
"suffienctly small" ds, but near bifurcation points only in cones around the branch. Thus, near 
a bifurcation point it is often not useful to choose very small ds. To circumvent this and similar 
problems we provide the function pmcont. The basic idea is explained in Fig. 11. 

Instead of using just one predictor (u^, X^) = (uq, Aq) + p.dsr, pmcont creates in every contin- 
uation step the predictors 

{u^, A^) = {uq, Aq) + i p.ds T, z = 1, . . . , p.mst, 

and starts a Newton loop for each, pmcont then monitors the convergence behaviour of each loop 
to decide whether it yields a "good" point, i.e., a point on the present branch. The criterion is that 
in each Newton step the residual has to decrease by a factor < a < 1, i.e., 

||GK+i,A„+i)|| < a||GK,A„)||, (35) 

Still, the bifurcation diagram in Fig. 10 is computationally quite expensive (about 40 minutes on a quad-core 
desktop PC, with about 60.000 triangles on average), and therefore the init-function p=schnakinit (p,m,n,nx,del) 
takes the domain sizes m, n, the deformation parameter 6 and the startup spatial discretization nx as parameters, 
schnaklldemo .m (or schnakllcmds.m) then uses m = n = 1 and 6 = 0.97 and only takes a few minutes for a 
bifurcation diagram similar to Fig. 10 over the smaller domain, while schnak22cmds .m treats the m — n — 2 domain 
in Fig. 10. 
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Algorithm pmcont 

1. Multipredictors. {u\ A*) = {uq, Aq) + i p.ds r, i = 1, . . . , p.mst 

2. Newton-loops (parallel). Use (35) to identify "good points". 

3. Tangents (sequentially). Calculate new tangents ri, . . . ,rm at good 
points, using (12). 

4. Bifurcation detection and localization (parallel). 

5. Postprocess (sequentially). Call ufu, save, plot, return to 1. 

Figure 11: Sketch of the basic idea of multiple predictors and convergence monitoring, and pseudocode 

(2) 

of pmcont. The arrows from, e.g. "2ds" to just illustrate the result of the Newton loops, not the 
hyperplanes {{u, A) G M^"'p+^ : (tq, {u, A))^ = ds} as in Fig.l. 

otherwise the loop is stopped. The heuristic idea is that if the Newton loop converges slowly, then 
probably the solution is on a different branch, because the loop has to (slowly) change the solution 
"shape". 

For the crucial parameter a, which describes the desired convergence speed, we recommend 
trying a = 0.1. Of course, these heuristics in no way guarantee that no branch switching occurs, 
or anyway that we get convergence for long predictors (ix^ A^) = (uq^Xq) + i • dsr with i > 1. 
But in practice we find the idea to work remarkably well. See also §5.2 for an example of the 
"unreasonable effectiveness" of pmcont, together with an example that long predictors in pmcont 
tend to branch-switching in imperfect bifurcations. 

Concerning data structures, we need three additional parameters (with p the problem struc- 
ture): the number of predictors p.mst, a =p.resfac, and p.pmimax. These are used to gain some 
flexibility for (35) via stepsize control. If p.mst equals the number of solutions found and p.ds is 
smaller than p.dsmax/p.dsincf ac, then p.ds will be increased by the factor p.dsincfac, essen- 
tially as in cent. If no solution is found and p.ds is greater than (l+p.mst)-p.dsmin, then the 
step size p.ds will be divided by 1+p.mst in the next continuation step. Finally, if p.ds is less 
than (l+p .mst) -p . dsmin and p.pmimax is less than p.imax, then p.pmimax will be increased to 
p.pmimax+1, and ||G('^^+p.pmimax, A^+p.pmimax) || < p.resfac||G('^^, A^)|| is required instead of (35). 
The only new functions used in pmcont are pmnewtonloop .m and pmbif detec .m. 

Another advantage of the p.mst predictors is that the Newton loops can be calculated in 
parallel, which on suitable machines gives substantial speedups.^^ All final Newton iterates with 
a residual smaller than p. res are taken as solutions, and plotted and saved as in cent. Next, the 
tangents ri, . . . , r^n are calculated sequentially, because t/+i needs r/, and afterwards the bifurcation 
detection and localization is again in parallel. The last solution and its tangent will be used for the 
next continuation step. Mesh adaption/refinement is inquired at the start of pmcont, i.e., before 
generating the predictors, but not on the individual correctors. 

In summary, for p.resfac=l and p.mst=l we have that pmcont is roughly equivalent to cont, 
except for slightly less versatile mesh adaption and error estimates. For p.mst> 1, pmcont takes 
advantage of parallel computing, and is often useful to avoid convergence problems and undesired 
branch switching close to bifurcation points. The main reasons why we (currently) keep the two 
version and do not combine them into one is that cont is simpler to hack and implements Keller's 
basic, well tested algorithm. 

^^However, for instance on the "hot hexagon branch" h in Fig. 10 we need to use a = 10~^ to avoid branch- 
switching, see schnak22cmds .m. 

^^Here we use the Matlab Parallel Computing Toolbox in an elementary setup with basic monitoring of open 
kernel threads. 
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5 Three classical examples from physics 



In this section we consider models for Bose-Einstein (vector) solitons, Rayleigh-Benard convection, 
and the von Karman system for buckhng of an elastic plate, as examples for systems with more 
than two components, and with EC implemented via gnbc as described in §3.1.4. The largest and 
most complicated system (in the sense of number of components and implementation of EC) here 
is the von Karman system. Hence, for this we also explain the coding in pde2path in most detail, 
while for Bose-Einstein solitons and Rayleigh-Benard convection we mostly refer to the m-files for 
comments. 

5.1 Bose-Einstein (vector) solitons (gpsol) 

As an example with x, y dependent coefficients, nontrivial advection, and interesting localized 
solutions we consider (systems of) Gross-Pitaevskii (GP) equations with a parabolic potential that 
arise for instance as amplitude equations in Bose-Einstein condensates. 

5.1.1 The scalar case 

First, following [18] we consider the scalar equation 

idt^l^ = + - ^IV^I V, (36) 

where = ip{x^y^t) E r'^ = x'^ + y^, and a = 1 (focussing case). This has a huge number 
of families of localized solutions, aka solitons, which may be time periodic, standing or rotating 
in space. Going into a frame rotating with speed uj and splitting off harmonic oscillations with 
frequency /i, i.e., 

i^{x,y,t) = ^r,<l>-ut)e-'^'\ (37) 

we obtain 

[d^ + ^dr + - icjde + /i - r^] $ + a\^\^^ = 0. (38) 

A typical ansatz for (approximate) solutions has the form 

$(r,^) = A(l){r/a){cos{m9) + ipsin(m^)), G M, e.g. (/)(p) = p'^L^^\p^)e-P^/^, (39) 

with L^^^ the n^^ Laguerre polynomial. Plugging this into (38) yields expressions for A, a,^, 
for approximate solutions. The case n=0 and p=0 corresponds to so called (nonrotating, since 
a; = 0) real m-poles, and \p\ = 1 to a so called radially symmetric vortex of charge m, while 
the intermediate cases < |p| < 1 give to so called rotating azimuthons with interesting angular 
modulations of |$|. 

Our goal is to calculate these solutions numerically with pde2path. Returning to cartesian 
coordinates, i.e., setting $(r, ^) = u[x^y) + \v{x^y) we obtain the 2-component real elliptic system 

- /\u+ (r^ - ii)u - \U\^u - uj{xdyV - yd^v) = 0, (40a) 
-Av + (r^ - fi)v - \U\^v - uj{yda:U - xdyu) = 0, (40b) 

where = u^+v'^. Our strategy is to use (39) for a; = and to continue in A := cj. A measure for 
the deformation of multipoles into vortices for the numerical solutions is the "modulation depth" 
p of the soliton intensity 

p = max |Im$|/ max |Re$| = max \v\/ max \u\. (41) 
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The result of typical continuation of a quadrupole using a stiff-spring approximation of DBC for v 
on domain = [—5,5]^ is shown in Fig. 12, (a)-(f), see gpf .m, gpjac.m, gpcmds.m and gpinit.m, 
and also plotsol.m in directory gpsol for the customized of plotsol. 
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Figure 12: (a) Continuing a quadrupole to an azimuthon to a 2-vortex, p = max max |v| over 
A = cj, (p branch, black); and continuing a vector dipole to a vector-azimuthon to a vector-vortex, 
Pi = max I I'll /max 1 1^1 1 over X = uj {q branch, red), ja = 2 resp. /ii = 2,/i2 = 2.2. (b),(c) vectorfield 
plots at A = (quadrupole) resp. A ^ 0.66 (azimuthon). (d)-(f) \U\ and arg('u + i'^) as indicated. (g),(h) 
vectorfield plots for the first condensate t/ji as indicated, the second condensate is similar. Also see the 
customized plotsol.m in gp. The phase plot in (f) is somewhat ragged due to the coarse mesh away from 
the center. For the p branch we used fixed p.nt=5208 from mesh-refinement during init. In the q branch 
we used q.ainod=8 with q.inaxt=9000 which lead to q.nt=9500 at, e.g., point 10. 

The following remarks are in order. First, we switch off stability or bifurcation tracking 
(spcalcsw=0, bif checksw=0), see Remark 5.1 below. Second, the linearization of (40) is given 
by (recall that A = cj) 



-A 






-A 



+ 



-/X + - 3u^ — v'^ 
2uv 



-/i + r 



—2uv 



uv 



-A 







-xdy + ydx 



y yda 





Thus, the last term is a good example how to use assemadv with a relatively complicated b. Third, 
some problems should be expected from the large number of solutions of (38), in particular the 
phase-invariance: if $ is a solution, so is for any a, or equivalently, (40) is invariant under 

cos a — sin a 



multiplication with ( """"" """""" i Thus, even for all parameters fixed, solutions of (40) always 
^ Vsma cos a y ' ^ v / ^ 

come in continuous families, and thus Gu as a linear operator on [L^(M^)]^, say, always has a zero 
eigenvalue. See also [18] and the references therein for some tricks for the numerical solution of 
(38). Rather remarkably, we need none of these tricks, presumably since numerically in Gu the zero 
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eigenvalue is perturbed sufficiently far away from 0. One trick we do use is to start with a coarse 
mesh of 30 x 30 points, first take some (rather arbitrary) monopole as dummy-starting guess, use 
meshref to generate a rather fine mesh in the center, define a quadrupole initial guess using (39) 
on that first refined mesh, and then refine again, yielding the (still small) number of 5208 triangles 
for this continuation, with an error-estimate less than 0.01. On a small laptop computer the whole 
continuation takes about a minute. 

5.1.2 A two— component condensate 

The above can be generalized to multi-component condensates [19]. For two components, we then 
have coupled GP equations of the form, e.g., 

i^t^i = [-A + - - ^^i2|^2p]^i, i5t^2 = [-A + - (j\^2? - 92i\^i\^]^2, (42) 

where ^12,^21 are called interspecies interaction coefficients. Physically, it makes sense to use 
ansatze of the form (37) with different ji but equal cj, i.e., '0j (x, t) = ^j{r^ (p — ujt)e~^^^^. Next we 
can use the form (39) for each component $j and classify the thus obtained approximate solutions as 
soliton-soliton, soliton-vortex, soliton-azimuthon etc pairs. To calculate such solutions numerically 
we set $j(r, ^) = Uj{x,y) + Vj{x,y) and obtain an elliptic system of the form (40) but with four 
real equations. This has been implemented in vgpf , with Jacobian vgpjac. Figure 12 (a),(g),(h) 
shows the continuation of a two-dipole obtained from vgpcmds. Similar remarks as for the scalar 
case apply, see the comments in vgpcmds. m. 

Remark 5.1. This clearly was just a very introductory demo of continuation of solutions of (36) 
resp. (42); there are many more and interesting branches, and further questions, again see, e.g., 
[18, 19]. Interesting questions concern, e.g., the dependence of vector solitons on — /i2| which can 
for instance be studied by fixing oj and continuing in A = /i2, or the effect of including a periodic 
potential, leading to gap solitons [6]. Some of these questions will be considered elsewhere. J 

5.2 Rayleigh-Benard convection (rbconv) 

As an example from fluid dynamics we consider two-dimensional Rayleigh-Benard convection in 
the Boussinesq approximation in the domain Q = [—2,2] x [—0.5,0.5]. In the streamfunction 
formulation the stationary system reads 

-Aip + oj = 0, 

—aAu — crRdxO + dx^l^dzUJ — dzi/^dxco = 0, (43) 
- - dx^ + dx^d^e - dz^dxO = 0, 

with streamfunction 7/^, temperature ^, and the auxiliary u = A^^. Moreover, a is the Prandtl 
number, set to 1 here, and R the Rayleigh number, which will be the continuation parameter. The 
implementation of (43) in pde2path is relatively straightforward, including analytical Jacobians, 
see rbconvf .m and rbconv jac.m 

The boundary conditions at the top and bottom plates are taken at constant temperature and 
with zero tangential stress 

= dzz^ = 6> = 0, at z = ±0.5; 

note that this means Aip = at z = ±0.5. Motivated by the analysis in [14], laterally we consider 
on the one hand "no-slip" (and perfectly insulating) BC 

^ = = dxO = 0, at ;z: = ±L, (44) 
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and on the other hand "stress free" BC 



'0 = dxx^ = dxO = 0, at X = zbL. 



(45) 



See the comments in rbconvbc_noslip .m resp. rbconvbc_stressf ree .m for the implementation 
(approximation) of these BC based on (2). 

In both cases it is known that continuation of the trivial zero state for increasing R gives a 
sequence of bifurcations alternating between even and odd modes. The stability thresholds are 
plotted in [14], Figure 1(a) for (45) and 1(b) for (44). We use these to choose initial values of 
A = i? for the first two bifurcations, respectively. 

For the no-slip case (44), the resulting bifurcation diagram is plotted in Fig. 13, which corre- 
sponds to the sketch Fig. 2 in [14]. No secondary bifurcations are found up to i? = 900. 



(a) Bifurcation diagram (b) qlO (top) and r8 (bottom) 

3r 




\ \ \ 



Figure 13: (a) Partial bifurcation diagram of (43) with (44). (b) sample solutions {ip^ and arrows indicating 
the fluid flow) from (a). See also arrowplot .m for how to produce the quiver plots. 

For stress- free BC we obtain the bifurcation diagram in Fig. 14(a), which corresponds to the 
case b^^ > a'^ ^h' > in Fig. 3 of [14]. Here the secondary symmetry breaking pitchfork from [14] is 
turned into an imperfect pitchfork. The x —x reflection symmetry is broken by the triangle data 
of the mesh (here we use poimesh) and the stiff- spring approximation of the boundary conditions. 
We have located the stable branch s of the imperfect pitchfork by time-integrating with tint 
from the unstable branch in the suitably chosen unstable direction. 



(a) bifurcation diagram 

3r 



(b) sfl-40 (top) and sf2-30 (bottom) (c) sf2-60 (top) and sf3-65 (bottom) 
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Figure 14: (a) Partial bifurcation diagram of (43) with stress-free b.c. and sample solutions {ip) from 
(a). Here sfl-40 and sf2-30 are approximately symmetric, while sf2-60 and sf3-65, generated in a 
(numerically) imperfect pitchfork around R = 860, are not. 

The demos run on a rather coarse mesh of 100 x 25 gridpoints, because even with assembled 
Jacobians the calculations are rather slow due to a non-simple structure of the Jacobians. Thus, 
we use pmcont for the bifurcating q and r branches, which gives a huge speed advantage and works 
remarkably well even directly after bifurcation from the trivial branch, where long predictors are 

^^which is not equivalent to the time integration of the time-dependent Boussinesq equations 
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far off the actual branches. ^ In Fig. 14, however, we switch back to cont when approaching the 
imperfect pitchfork since pmcont tends to switch from the r to the s branch via long predictors. 

Increasing the number of meshpoints brings the diagram closer to a symmetry breaking pitch- 
fork. In fact, for both boundary conditions, qualitatively the same bifurcation diagram can be 
found for rather coarse meshes, but the location of the branches can be off by order 100 in R. 

5.3 Von Karman description of the buckling of plates (vkplate) 

The von Karman equations 

-A^v - Xdlv + [v, w] = 0, -A^w - ^ [v, v] = 0, (46) 

can be derived to describe the deformation of an elastic (rectangular) plate = [—/a;, Z^;] x [— /y, /y] C 
under compression. Here : -> M is the out of plane deformation, t^; : J7 — > M is the Airy 
stress function, = [d^ + 5^)^ is the squared Laplacian, A is the compression parameter, and the 
bilinear form [•, •] is given by 

There are a number of choices for the boundary conditions for (46). For v one can choose for 
instance between (in the notation from [10]) 

I(v) : V = Av = on (simply supported), 

II(v) : V — Av = on y = ±/y, v — dnV = on x = zb/^;, 

(simply supported on the sides, clamped at the ends) 
III(v) : V = dnV = on (clamped on whole boundary). 

Similarly, for w we may consider, on dfl^ 

I(w) : w = Aw = 0, II(w) : dnW = dn{Av) = 0, III(w) : w = dnW = 0. 

Clearly, for all BC-combinations and all A the trivial state = it; = is a solution. Mathe- 
matically, the combination a) (I(v),I(w)) (sometimes as a whole called simply supported) is most 
simple because it allows an easy explicit calculation of bifurcation points from the trivial branch. 
However, [27] argues that physically the combinations b) (II(v),I(w)) or c) (II(v),II(w)) are more 
reasonable, and various combinations and modifications have been studied since, see [4] and the 
references therein for an overview. 

Here we focus on case b) since this yields secondary bifurcations, called "mode jumping" in 
this field. The other cases can be handled quite similarly and, e.g., a) is in fact slightly simpler. 
The aim is to show how (46) can be put into pde2path and thus recover a number of interesting 
bifurcations. 

Clearly, the first idea to set up (46) would be to introduce auxiliary variables Av^ Aw and set 

U (ui, U2, Us, U4) = {v, Av, w, Aw) 

to obtain the (quasilinear elliptic) system 

I -A 1 \ 

-Xdl -A 

-A 1 

V -a) 

^^E.g., in Fig. 13 all points are calculated from the essentially "vertical" predictors at bifurcation! To check that 
this does not miss any (possibly imperfect) bifurcation we compared with cont with small ds and obtained the same 
branches but much slower. 



/ \ 



\\[uuui\) 



0, 
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for instance in case a) with homogeneous Dirichlet BC ui — U2 = us = U4 = 0. The problem with 
this formulation in pde2path are the derivatives d^ui, . . . , dxdyU^ in the nonlinearity. In principle, 
these can be obtained from calling pdegrad, pdeprtni, and pdegrad again^^. However, the first 
problem is that this introduces some averaging into the second derivatives, in particular at the 
boundaries. The second problem is that with this approach we have no easy way to generate 
the Jacobian of G since pdegrad/ pdeprtni neithers fit to matrix assembling nor to numerical 
differ ent iat ion . 

Thus, here we choose to introduce additional auxiliary variables, i.e., set 

u = (i;, Ai;, Aw, d^v, dyV, dxdyV, d^w, dyW, dxdyw) G R^^. 

For instance, = d^ui can then be simply added as a linear equation —d^ui + 1^5 = in the 
pde2path formulation^^. However, since this way we get a number of indefinite equations, in 
particular the mixed derivatives —dxdyUi + 1^7 = and —dxdyUs + uio = 0, here we use an ad hoc 
regularization and set —d'^ui + (1 — 5A)u^ = with small 5 > (i.e., 6 = 0.05 numerically) and 
similarly for ne, . . . , i^io- Thus, instead of (46) we now really treat the problem 

-A'^v - Xdlv + [SvxxSwyy - 2SvxySwxy + SvyySwxx) = 0, ^^^^ 

-A'^W - {SVxxSVyy - SVxySWxy) = 

with the smoothing operator 5 = (1 — 5A)~^ . However, for small 5, comparison of our results with 
the literature shows that the regularization plays no qualitative or even quantitative role (in the 
parameter regimes we consider). 

Thus, we now have a 10 component system, and to illustrate its implementation in pde2path 
we write it in the form [—C + A)u — f = Q with 

/ = (0, —{u^u^ — 2uyUio + uqUs)^ 0, u^uq — u^, 0, 0, 0, 0, 0, 0)"^, and 



-C + A = 



( -A, 



%2 



-dy21 



V 



— 0x109 

— dyll3 

-dxdyiij 



V 



Here, (for layout reasons) V — bA + 1, and the subscripts 1, 5, 17, . . . denote the starting positions 
of the respective 2x2 tensor stored in the "400 rows vector" c.^^ The superscripts 11,33,... 
denote the positions in the "100 rows vector" a, and for T) subscripts refer to 6 A and superscripts 
to +1. See vkf .m. Similarly, it is now rather easy to put the linearization into pde2path, i.e., 
the second and fourth row of as a 10 x 10 matrix read 



2nd row: 
fu, 4th row: 



(0 
(0 



^6 



54 



52 



62 



2^/10" 



-Uq 





72 



-^5 




82 







^e.g. , [ulxt , ulyt] =pdegrad (p . points , p . tria , u ( 1 : p . np) ) ; ulx=pdeprtrLi (p . points , p . tria , ulxt ) ; 

[ulxx,ulxy]= pdegrad (p. points, p. tria,ulx) ; could be used to calculate (approximate) dxUi 

^°For the latter the next-next-neighbor effect of pdegrad/ pdeprtni does not comply with the Jacobian stencil, 
as explained in Remark 3.1. 

^^see below for the BC for 1^5, ... , i^io 

^^I.e., — Ai yields ci = [1; 0; 0; 1] stored in positions 1 to 4 in c, dx^ yields C2 = [1; 0; 0; 0] stored in positions 5 to 
8 in c, and so on. 
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Here again the superscripts give the positions in f u. Of course, the full fu — fu — a also contains 
the constant coefficient terms at positions 11, 33, 45 etc from A] see vkjac.m. 

It remains to encode the boundary conditions. First note that v — and w — Aw = imply 



'^5 — "^xx - 
U8 



and uq = Vyy = on horizontal edges, and 
on vertical edges, but no condition for - 
- and uq = Wyy = on all edges. 



and 



For on the vertical edges and u^^uiq on all edges we take homogeneous Neumann boundary 
conditions. To put this into pde2path via (2) we thus need two boundary matrices and g^. For 
the horizontal boundaries {y = i/y), has diagonal q^^i^ssssssOssO). For the 
vertical boundaries {x = ±lx) q^ has diagonal q'^^^OOssOsOssO) and additionally 
^2 1 = 5, where s = 10^ stands for the stiff spring constant. Positions 7 and 10 in q^ and q^ give 

0^ 



the Neumann BC for u^^uiq, while the top left 2x2 block ( ^ ^ ) in gives dnUi 



via the 



first row and U2 = via the second row. 

The (analytical) calculation of bifurcation points from (i;, i^;) = in case b) is rather tedious, see 
[27]. There, motivated by mode-jumping, the particular interest is in (the lowest) double bifurcation 
points, which yields / = \/k{k + 2) with eigenfunctions wi{x, y) = (^^sin(fcj) — sin((/c + 2) j)) sin(^) 
resp. W2{x,y) = (cos(A:j) — cos((A: + 2) j)) sin(y) (over the domain [0,/7r] x [0,7r]). The first bifur- 
cation is then obtained for A: = 1, hence / = \/3. The idea is to perturb / slightly which may lead 
to secondary bifurcations between branches coming originally from the same A. 

Putting all these ideas together we indeed get a secondary bifurcation between the first two 
primary branches, see Fig. 15. A number of further bifurcations from the trivial branch is also 
detected and can be followed. However, in the tutorial run vkcmds we use a rather coarse mesh 
with 1250 triangles, which should be refined before following higher bifurcations. 





qbpi^ 


1 













-1 



at rbpl 



u at rbpl 



at wp3 



at wp6 



I 





Figure 15: Secondary ("mode jumping") bifurcation (w-branch, red) in the (regularized) partially clamped 
plate (47): maxiii over A, selected plots of ui and us from the bifurcation diagram, and ri at the third 
bifurcation (A ~ 9.1) from the trivial branch, ly = 7r/2, = 47r/5, regular mesh with 25 x 25 points (1250 
triangles). Error estimate 0.3 at, e.g., rbpl. By mesh refinement we can obtain an error-estimate~ 0.045 
with nt= 14916. Then, however a typical step takes a couple of minutes, where about 80% of the time is 
spent in blss or Iss (standard setting). We expect that this can be optimized considerably, but here we 
content ourselves with the "proof of principle" setup for the 10 components system for (47). 
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6 Discussion 



Clearly, numerical continuation and bifurcation analysis for 2D elliptic systems poses additional 
challenges compared to algebraic equations or ID BVP, partly of course due to the more demanding 
numerics, but in particular also due to the typically very rich solution and bifurcation structure. 
With pde2path we believe to provide a general tool that works essentially out-of-t he-box also for 
non-expert users and allows to start exploring such systems and the rich zoo of their solutions. Of 
course, in many respects this is just a first step, and probably the main entries on our to— do— list 
are: 

1. Implement some more general bifurcation handling, in particular bifurcation via nonsimple 
eigenvalues as these are quite ubiquious in 2D systems due to various symmetries. 

2. Implement some (genuine) multi-parameter continuation. For instance, the bifurcation to 
travelling waves generically requires a second parameter 7 (the wave speed) to adapt, and 
consequently we need to further extend the "extended system" (8) by one more equation, 
the "phase condition" .^^ We believe that our set-up of pde2path is sufficiently modular and 
transparent such that this and similar adaptions will pose no implementation problems, but 
for now we confine ourselves to the basic one-parameter continuation and simple bifurcations. 
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